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Preface 


Nuclear thermal hydraulics is the application of thermofluid mechanics within the nuclear indus- 
try. Thermal hydraulic analysis is an important tool in addressing the global challenge to reduce 
the cost of advanced nuclear technologies. An improved predictive capability and understanding 
supports the development, optimisation and safety substantiation of nuclear power plants. 


This document is part of Nuclear Heat Transfer and Passive Cooling: Technical Volumes and Case 
Studies, a set of six technical volumes and four case studies providing information and guidance 
on aspects of nuclear thermal hydraulic analysis. This document set has been delivered by Frazer- 
Nash Consultancy, with support from a number of academic and industrial partners, as part of 
the UK Government Nuclear Innovation Programme: Advanced Reactor Design, funded by the 
Department for Business, Energy and Industrial Strategy (BEIS). 


Each technical volume outlines the technical challenges, latest analysis methods and future direc- 
tion for a specific area of nuclear thermal hydraulics. The case studies illustrate the use of a subset 
of these methods in representative nuclear industry examples. The document set is designed for 
technical users with some prior knowledge of thermofluid mechanics, who wish to know more about 
nuclear thermal hydraulics. 


The work promotes a consistent methodology for thermal hydraulic analysis of single-phase heat 
transfer and passive cooling, to inform the link between academic research and end-user needs, 
and to provide a high-quality, peer-reviewed document set suitable for use across the nuclear 
industry. 


The document set is not intended to be exhaustive or provide a set of standard engineering ‘guide- 
lines’ and it is strongly recommended that nuclear thermal hydraulic analyses are undertaken by 
Suitably Qualified and Experienced Personnel. 


The first edition of this document set has been authored by Frazer-Nash Consultancy, with the 
support of the individuals and organisations noted in each. Please acknowledge these documents 
in any work where they are used: 


Frazer-Nash Consultancy (2021) Nuclear Heat Transfer and Passive Cooling, 
Volume 3: Natural Convection and Passive Cooling. 
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1.1 


Introduction 


This technical volume considers fluid flow and heat transfer in passive cooling systems, focusing 
primarily on flows driven by buoyancy. Such flows are of particular relevance to passive cooling 
systems in the nuclear industry and present particular technical challenges. 


This volume is part of a set of technical volumes, and it is recommended that Volume 1 (Introduction 
to the Technical Volumes and Case Studies) and Volume 2 (Convection, Radiation and Conjugate 
Heat Transfer) are reviewed for context. This introduction defines and contrasts the characteris- 
tics of passive cooling systems with active cooling systems, introduces natural, forced and mixed 
convection and presents typical applications of passive cooling systems. 


Passive and Active Cooling Systems 


IAEA (1991) defines a passive safety system as: Either a system which is composed entirely of 
passive components and structures or a system which uses active components in a very limited 
way to initiate subsequent passive operation. Therefore, passive safety systems are designed to 
reduce or remove the need for active intervention by an operator or control system to bring and 
maintain a reactor to a safe shut-down state, in the event of a particular scenario occurring. 


Passive cooling systems therefore transfer heat without needing external inputs (at least after sys- 
tem operation is initiated), and typically take advantage of natural forces or phenomena such as 
gravity, buoyancy, pressure differences, conduction, thermal radiation or natural heat convection to 
accomplish safety functions without requiring an active power source (IAEA, 2009a). 


Active Cooling Systems: Most of the cooling systems in current Nuclear Power Plants (NPPs) 
are closed loop active systems, and they are often numerous, partly because indirect cool- 
ing using chains of systems connected to each other is common. These systems, their 
support systems! and all the associated components must be designed and integrated 
within the whole NPP, manufactured in a supply chain, installed, commissioned and main- 
tained through life. As such, while designing and testing these systems may appear relatively 
straightforward from a technical perspective, they may be expensive and present risks during 
construction and operational life. 


Passive Cooling Systems: By contrast, passive cooling systems do not contain components 
such as pumps and fans, and may use simpler cooling chains or no cooling chains at 
all. Therefore, new reactor designs are making more extensive use of passive safety fea- 
tures because they are intended to achieve the same or higher reliability using fewer sys- 
tems/components (less complexity), thus reducing capital and Operation and Maintenance 


1 Such as electrical, control and instrumentation, motor cooling, other cooling systems and all their support systems. 
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(O&M) costs associated with the installation and maintenance of mechanical, control and in- 
strumentation support systems. For example, a reactor design developed by NuScale* uses 
natural circulation of the reactor coolant flow around the primary circuit (e.g. without main 
coolant pumps) under normal operation and loss of power conditions. Examples of passive 
cooling systems include core cooling, containment cooling, safety injection and decay heat 
removal systems. 


Overall, IAEA (2002b) notes that using passive safety systems is a desirable method of increasing 
simplification and potentially reliability, and that passive systems should be used where appropri- 
ate. However, the nature of the flows, complex coupling between the flow and temperature field, 
and weak or unstable driving forces associated with passive cooling systems (e.g. natural circula- 
tion) can make justification of the plant’s operation and safety performance more challenging: 


« Can result in physically larger system components, because the flow velocities associated 
with passive cooling systems are normally lower than active cooling systems. 


¢ Might require more detailed analytical modelling, experimental work and testing to design, 
substantiate performance and confirm reliability, particularly where a graded approach indi- 
cates a high level of scrutiny is appropriate (Volume 1, Section 2.2.5), because: 
— The cooling flow rate and heat transfer performance are coupled (although a benefit of 
this is that the cooling flow rate may increase as heat is transferred to the flow). 


— Actuation may not be as instantaneous as switching on an active system, so more time 
may be needed to reach full cooling performance. 


— Thermal mass may be more significant than would be the case for active systems. 


— Natural or mixed convection flow fields are likely to occur, and these are often spatially 
complex and intrinsically unsteady. 


— System flows are more likely to be at transitional Reynolds numbers, due to the lower 
driving forces associated with natural circulation. 


IAEA (1991) also recognises that passive cooling systems may require modest control functions 
to enable operation. Some guidance on categorising the ‘passivity’ of passive systems from Cat- 
egory A (more passive) to Category D (less passive) is provided, and discussed further in IAEA 
(2002b). This considers aspects such as the level of control functionality needed and the impor- 
tance of moving parts or fluids. For example, a static wall used to segregate systems from one 
another might be Category A, a system with moving fluids but no moving parts, control functions 
or support systems might be Category B, and a system with active initiation followed by passive 
execution (e.g. that opens valves to connect a natural circulation loop, see Section 1.2) might be 
Category D. This illustrates that the overall robustness of a passive cooling system may not begin 
and end with thermal hydraulics. The wider system must be appropriately designed as a whole to 
provide the overall robustness or reliability required. 


As noted above, achieving the potential benefits of passive cooling may involve more sophisticated 
technical design work than might be needed for active cooling systems. A motivation for this tech- 
nical volume is to help develop understanding of how sophisticated this work might need to be and 
how it might be approached. 


2 www.nuscalepower.com 
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1.2 Natural, Forced and Mixed Convection and Circulation 


Convection: The terms natural, forced and mixed convection describe mechanisms of convec- 
tive heat transfer, rather than overall system behaviour: 


* Natural convection?: Local heat transfer arising from flow driven by local buoyancy effects. 
* Forced convection: Local heat transfer arising from flow driven by work input. 


¢ Mixed convection: Local heat transfer arising from aspects of both natural convection and 
forced convection. 


Flows that may result from natural, forced and mixed convection are illustrated in Figure 1.1. 


MAA # #£3Moory WATATATATATAVAT 


(a) Natural Convection (b) Forced Convection (c) Mixed Convection 
Figure 1.1: Typical flows associated with natural, forced and mixed convection. 
Circulation: The terms natural, forced and mixed circulation describe what is driving the flow 

from the perspective of the overall system: 


¢ Natural circulation: Flow in the system as a whole driven by buoyancy. 


* Forced circulation: Flow in the system driven by mechanical equipment like a pump or fan, 
as seen in most active cooling systems. 


¢ Mixed circulation: Flow in the system driven by both buoyancy and mechanical equipment 
(which may be aiding or opposing each other). 


Illustrations of natural, forced and mixed circulation are shown in Figure 1.2. 


Binnie 


(a) Natural Circulation (b) Forced Circulation (c) Mixed Circulation 


Figure 1.2: Example systems containing natural, forced and mixed circulation. 


3 The term ‘free convection’ is sometimes used to describe natural convection. The term ‘free’ normally refers to flows that 
are away from walls, and as these can inhibit natural convection this term may be appropriate (e.g. in pools). However, to 
avoid confusion with other ‘free’ flow features, the term ‘free convection’ is not used in this volume for clarity. 
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For a system operating using forced circulation, such as a pumped primary circuit or auxiliary 
cooling system under normal operation or faults where the pumps are running, most of the flow 
fields in the circuit, pipework, valves and pumps are likely to be forced convection. However, mixed 
convection could be encountered if there are locations where relatively low forced flow velocities 
occur near strong thermal gradients. In addition, natural convection may occur in areas where the 
forced flow is quiescent (such as inside ‘dead’, or isolated, pipework spurs). 


Mixed circulation is unusual, but could be encountered where both mechanical components and 
temperature gradients are significant in driving the flow. For example, a heavy pump slowly decel- 
erating to a stop in a system containing large temperature gradients (as might occur in a primary 
circuit if the main coolant pump electrical supply is interrupted). 


For natural circulation, natural convection is likely to be encountered where there are strong thermal 
gradients (such as at the start of a transient or in a pool, where flow is likely to be quiescent). 
However, where natural circulation is well established in a piped system, the flow in the vicinity of 
these thermal gradients may more closely resemble mixed convection as a result of the momentum 
of the flow in the system, even though this flow is entirely driven by buoyancy effects. The flow in 
parts of the system where thermal gradients are low (such as pipework) may resemble flow fields 
associated with forced convection, but unlike in forced circulation the overall flow rate and heat 
transfer performance are coupled, and flow velocities are likely to be lower. 


Passive cooling systems are likely to make use of natural circulation. However, as noted above, the 
resulting flow fields inside the system may resemble natural, mixed or forced convection. Natural, 
forced and mixed convection flows are considered in more detail in Section 2.2. 


Passive Cooling Applications 


Passive cooling systems are widely used in a range of applications. In NPPs, typical uses for 
passive cooling systems include: 


Cooling the primary circuit, particularly for decay heat removal in accident scenarios, to main- 
tain the integrity of the fuel and/or primary circuit (typically the first and second barriers that 
confine nuclear material within the power plant). 


* Cooling containment buildings internally or externally, to manage internal pressure in the 
event of a primary circuit breach and maintain the integrity of the building (which is typically 
the third barrier). 


Cooling spent nuclear fuel (again to preserve the integrity of confining barriers). 


Cooling the air in buildings to manage internal temperatures (supporting the continued oper- 
ation of systems and access by operators). 


Passive cooling systems must be carefully designed in a holistic manner to ensure that the overall 
performance of the NPP is appropriate, and generally adopt one of the following arrangements: 


1. Volumes of cold water that can be fed directly into the NPP’s cooling systems. Such systems 
are generally driven by gas pressure in accumulators or gravity feed from elevated tanks, 
and provide a finite amount of cooling (the cooling stops when there is no more water left). 
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2. Volumes of initially cold water that can be used as a heat sink (rather than being injected 
into the NPP’s cooling systems), again providing a finite amount of cooling. Heat is typically 
transferred into this volume from the plant using natural circulation in pipework, or thermal 
radiation in high-temperature reactors. 


3. Cooling chimneys that are used as a heat sink by a natural circulation loop that removes 
heat from the plant. This provides a long term source of cooling, providing the chimney and 
cooling loop are designed to remain effective in the prevailing weather conditions. 


4. Ducts or louvres that exhaust heated air from rooms within buildings and allow cool air to 
enter, transferring heat by natural circulation. 


Some Advanced Nuclear Technologies (ANTs) are designed to use the same cooling systems 
for normal operation and safety scenarios within their design basis (and to minimise changes in 
operating states) although additional systems may also be included (e.g. to manage design exten- 
sion conditions). Natural circulation is often proposed to cool the reactor vessel, reactor pits and 
the containment, either directly or indirectly through a closed cooling loop. Passive coolers may 
also be immersed into the primary pool and connected to a closed cooling loop. High tempera- 
ture reactors are likely to make more use of thermal radiation for cooling than is normal in Light 
Water Reactor (LWR) designs due to their higher operating temperatures (particularly in safety 
scenarios), and these designs often use ambient air as the ultimate heat sink for passive cooling 
systems. In designs where fuel is incorporated into the coolant, melting plugs may be used to drain 
the fuel into cooling chambers. These aspects are discussed in more detail in Volume 5 (Liquid 
Metal Thermal Hydraulics) and Volume 6 (Molten Salt Thermal Hydraulics). 


Details on the application of passive cooling systems to LWRs are provided in IAEA (2009a), 
which provides a high-level overview of the passive systems typically employed (the annexes de- 
scribe how around 20 designs, including integral designs, incorporate passive safety systems). 
Further discussion on the incorporation of passive cooling systems into LWRs is presented in 
IAEA (2002b). A number of papers considering various aspects of passive cooling system design 
are presented in IAEA (1996) and detailed case studies on three reactor designs are presented 
in IAEA (2013c). An example of the passive cooling systems in an operating NPP is shown in 
Figure 1.3. 


Passive ventilation systems are also used, particularly for cooling plant buildings and facilities for 
storing spent fuel. IAEA (2018) and Appendix 3 of IAEA (2005) contain a brief discussion on 
passive ventilation systems specifically. From a technical perspective, the physical phenomena 
relevant to Heating, Ventilation, and Air Conditioning (HVAC) systems are similar to those for me- 
chanical cooling systems (illustrated by the example in Section 2.2.2). However, HVAC systems 
are likely to have lower temperature gradients, different types of requirements (e.g. confinement 
of radioactivity within the plant using filters, managing humidity and pressure levels in rooms etc.) 
and different components (using ducts, dampers, fans and chillers instead of pipes, valves, pumps 
and heat exchangers). Open loop passive cooling using natural convection within an HVAC system 
is generally called natural ventilation. Further detail on aspects of HVAC design are included in 
Haines and Myers (2010) and ASHRAE (2017). 
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PCS GRAVITY DRAIN 
WATER TANK 
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Figure 1.3: An illustration of passive cooling features employed by the Westinghouse 
AP1000 Nuclear Power Plant. Further detail is available at 
www.westinghousenuclear.com/new- plants/ap1000-pwr/safety. 

This image is copyright of Westinghouse Electric Company LLC, is used with permission, 
and is not covered by the creative commons license defined in the legal statement for 
the present document. 
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2.1.1 


DP 


Technical Context 


This section introduces some of the key technical aspects of the buoyancy affected flows that may 
occur in passive cooling systems. The key flow phenomena and theory are presented, fluid material 
properties are discussed and key modelling challenges are considered. 


This section does not consider all the phenomena that might be relevant to the design of a whole 
system’, but guidance in this area is developed through a series of three reports developed by an 
International Atomic Energy Agency (IAEA) Coordinated Research Project (CRP), (IAEA, 2005, 
2009a and 2012), which in turn build on a number of Organisation for Economic Co-operation 
and Development (OECD) Nuclear Energy Agency (NEA) reports on passive safety systems. Heat 
transfer by phase change? is also not considered in detail, as described in Volume 1. 


Flow Phenomena 


Key relevant flow phenomena such as plumes, stratification and surface effects are introduced 
below. In many practical engineering applications, different flow phenomena may interact with each 
other within a system® or with external ambient conditions (especially for open loop systems). 
These interactions may be simple and intuitive, but may also be very complex and their behaviour 
counter-intuitive even to skilled engineers. Analysis may be needed to identify and understand 
such interactions and may also enable them to be avoided or simplified by design if they have an 
undesirable impact (Section 3). 


Plumes 


When the density of a volume of fluid is reduced relative to fluid surrounding it, buoyancy forces 
act to push the less dense volume of fluid upwards. This results in the less dense fluid accelerating 
upward in a buoyant plume and more dense fluid being drawn in beneath to replace it. In natural 
convection, this reduced density is generally caused by increased fluid temperature, which may 
result from contact with a heated surface or heated fluid entering a cold volume from elsewhere. 


Figure 2.1 shows illustrative diagrams of free and bounded plumes, which are considered sepa- 
rately below, and Figure 2.2 shows these plumes in a flow loop. In this volume, plumes are generally 
described as if they are of low density (e.g. hot); high density (e.g. cold) plumes generally show 
the same behaviours in the opposite sense. Plumes caused by fluid species with different densities 
may show similar behaviours to the thermal plumes discussed below (Woods, 2010). 

Such as might be considered for Phenomena Identification and Ranking Table (PIRT) analysis, described in Volume 1. 
Such as evaporation from free surfaces, or boiling and condensation in pools, heat exchangers, steam generators, boilers 
or heat pipes (Rohsenow et al/., 1998; Reay et al., 2014). 


Examples include interactions between plumes (IAEA, 2005), between plumes and stratification, and between plumes, 
stratification and wall boundary layers. 
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Stratification 


Velocity Profile 


| 


(a) Free Plume (b) Wall Bounded Plume 


Figure 2.1: Illustrations showing free and wall bounded plumes with stratification. 


Free Plumes (Pools and Plena): For a ‘free plume’ rising through fluid away from walls and 
surrounded by a more dense fluid (e.g. in pools or large plena), the core of the plume will continue to 
accelerate until the buoyancy forces acting on it are counteracted by energy dissipation associated 
with its motion (Section 2.2.2). As the plume rises, this mixing with the cooler surrounding fluid 
will cause entrainment of surrounding flow into the plume, resulting in an increasing mass of rising 
fluid, but a cooling and hence decelerating plume core. 


If there is a large height of cool surrounding fluid, the plume will eventually dissipate (or ‘mix out’) 
into the surrounding fluid, losing all of its momentum and transferring all of its excess thermal 
energy to the surrounding fluid. If the plume encounters other plumes, flow features, walls or free 
surfaces, it will also interact with these, resulting in complex flow features and further mixing. If the 
temperature of the surrounding fluid rises, the buoyancy forces driving the plume will reduce and 
a vertical plume will therefore tend to decelerate and widen. Plume behaviour is discussed further 
in Gebhart et al. (1988), Rodi (1982) and Turner (1973). Rising temperature with height is often a 
result of stratification; as well as modifying plume behaviour, stratification can also decouple flows 
above and below an abrupt temperature gradient (Section 2.1.2). 


Where jets from pipework discharge into a large volume, the initial flow is likely to be driven by 
momentum effects. However, the momentum of the jet will decrease through mixing as it moves 
through the volume. If the density of the jet is different to the surrounding fluid (e.g. because it 
is hot), buoyancy effects are likely to become dominant as the jet develops and, given enough 
distance, the jet is likely to develop into a free plume (Gebhart et a/., 1988). 


Bounded Plumes (Loops and Channels): For a ‘bounded plume’ rising through fluid near ad- 
jacent walls within a system (e.g. in pipes within loops or channels within a reactor core), the initial 
flow phenomena are likely to be similar to a free plume with the addition of dissipation caused 
by wall-bounded flow phenomena like momentum and thermal boundary layers (Section 2.2.2). 
However, once the plume displaces a significant amount of fluid, wider system effects become im- 
portant. If the wider system is slow to react’, the flow adjacent to the plume may flow strongly in 
the opposite direction to the plume as a result of mass continuity, causing local mixing and slowing 
the development of the plume itself. Once plumes become able to rise up through the system, the 


This is likely, because the low driving forces in passive systems tends to lead to large reactors and components being used 
to reduce power density and flow pressure losses, both of which increase the mass of fluid in the system 
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Volume 3 


forcing effect of the heated flow on the system as a whole is likely to depend on the instantaneous 
(i.e. time-varying) balance between: 


1. What density variation is generated over what height (i.e. the magnitude of the overall hydro- 
static pressures driving the flow). 


2. The momentum of the heated flow (i.e. the mass of flow moving and its velocity, resulting 
from the acceleration due to the driving force). 
3. The energy dissipation due to viscous effects around the whole flow. 


4. The flow behaviour in the wider system, which may or may not be strongly coupled. 


An example of free and bounded plumes occurring in a flow loop is shown in Figure 2.2. In a 
closed loop system containing a heater and cooler, complex and coupled behaviour may occur at 
every length scale, from detailed local flow interactions to whole-system flow oscillations (or even 
reversal, see IAEA, 2012 and Section 2.4). This flow behaviour may be coupled with additional 
phenomena, such as geometry distortions or phase changes. 


Time 1 Time 2 


Flow Direction * 


Cooler 


Dimensionless 

Temperature 
Heater —_ Difference 

0 5 10 15 20 


Bounded plume 


10 15 20 


Free plume 


Figure 2.2: Free and bounded plumes occurring in a simple (two-dimensional, laminar) 
flow loop containing heated and cooled walls. 


Zele Stratification 
Stratified flows are often described as stable or unstable (Figure 2.3): 


« In stable stratification, a warmer (less dense) fluid sits stably above cooler (more dense) fluid. 
In this situation, slow molecular conduction is often the only mechanism to transfer heat from 
the warmer fluid into the cooler fluid below (in the absence of thermal radiation effects). As 
a result, heated fluid may collect under the free surface, sitting above the cooler fluid in the 
pool and preventing the fluid from circulating. 


* In unstable stratification, a cooler fluid is above a warmer fluid, leading to the development of 
buoyancy-driven flows. An example is Rayleigh-Bénard cells, a particular flow phenomenon 
with plumes arranged in a cell-like structure. 
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Rayleigh-Bénard 


Cells 
Cold 
(a) Stable (b) Unstable 


Figure 2.3: Stable and unstable stratification. 


Stratification is a particular feature of pools or volumes of fluid within systems, and can also be 
important in pipework. Some examples of scenarios where stratification in volumes or pipes may 
influence NPP design are illustrated in Figure 2.4: 


1. In systems where a cold tank of water is used as a heat sink, stratification can cause the 
capacity of this heat sink to reduce. This is because the tank effectively ‘fills down’ with 
warm fluid from the top, so cooling is likely to reduce once the interface between warm and 
cool fluid descends to the level of the heat exchanger transferring heat into the pool (IAEA, 
2013c). 


2. Stratification can cause sharp temperature gradients in large volumes (such as vessels in a 
primary circuit). lf this interface moves around (e.g. due to flow unsteadiness caused by natu- 
ral convection), components or walls that cross this interface will suffer repeated temperature 
fluctuations (thermal striping) that may cause thermal fatigue (IAEA, 2014 and Section 2.1.3). 


3. Stratification within system components can also cause the development of sharp tempera- 
ture fluctuations in the downstream flow (referred to as thermal ‘fronts’ or ‘slugs’), which may 
cause thermal fatigue or shock in downstream components (IAEA, 2012 and Section 2.1.3). 


4. Where a cold fluid enters a pipe full of slow moving hot fluid via a T-junction (such as safety 
injection into a Pressurised Water Reactor (PWR) primary circuit under natural circulation), 
the cold fluid may not mix with the hot fluid, but form a layer at the bottom of the pipe. This 
forms two regions of fluid, which can move in different directions with different velocities. 
The cold fluid may flow some distance down the pipe and cause Pressurised Thermal Shock 
(PTS), where components (such as parts of the Reactor Pressure Vessel (RPV) in a PWR) 
are exposed to high thermal stress (IAEA, 2012, IAEA, 2005 and CSNI, 2015a). 


5. In a spur of pipework or a pipe connected to a cooler containing stagnant flow, the heat 
transfer from the pipe may cause natural circulation within the pipe itself, with outgoing hot 
fluid forming a layer above the returning cold fluid (sometimes called a ‘thermosyphon’). This 
situation can cause large bending moments in pipework, problems associated with conden- 
sation of steam, and transport of hot fluid over long distances (which can cause an ‘induced 
steam generator break’ in a PWR primary circuit in two-phase conditions, CSNI, 2015a). If 
the pipe contains thermocouples these may be immersed in different areas of a complex flow 
field, making the readings difficult to understand (particularly for complex pipe routes). 


6. Where a region of pipework contains a low-point, this can trap cold fluid that has entered the 
system. The cold fluid may then remain in this ‘cold trap’ and resist or stagnate the natural 
circulation of hot fluid in the system (IAEA, 2005, Annex 15). It is necessary to understand 
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situations in which a passive cooling system might become ineffective, in order to set the 
Limits and Conditions for Operations (LCOs) for the NPP. 


Temperature Temperature 
Fluctuations Fluctuations 


3 


planes Blockage to Hot Flow 


6 


Figure 2.4: Example phenomena associated with stratification in volumes and pipes. 


Where there is a large, abrupt and stable density gradient (or ‘stratification interface’) the flows 
on either side of this gradient may become effectively decoupled, so that plumes or other flow 
features on one side may have a negligible impact on the flow on the other side (as in the cold 
trap example above). Where the density gradient is more modest, a warm plume or other flow 
feature may still be able to rise through the gradient, but the driving buoyancy forces will reduce, 
typically causing marked changes in the development of the flow field (Section 2.1.1, Woodcock 
and Dzodzo, 2000). The Atwood number, At (Section 2.2.1) may be useful in characterising the 
stability of such stratified flows. 


Surface Effects 


Walls: Heat transfer between solids and adjacent fluids, and the potential for thermal stresses to 
occur within the solid, is considered in detail in Volume 2 (Convection, Radiation and Conjugate 
Heat Transfer). The following aspects are characteristic of buoyancy driven flows: 


* These flows are often unsteady and poorly mixed (i.e. containing regions of relatively hot and 
cold fluid adjacent to each other). If such a flow is adjacent to a component or wall, this may 
suffer repeated temperature fluctuations that may cause thermal fatigue (often referred to as 
‘thermal striping’). 


* As noted in Section 2.1.2, stratification interfaces (whether stable or unstable) can cause 
temperature fluctuations that may cause thermal fatigue, or changes in temperatures within 
components that may cause high stresses. 


¢ Buoyancy effects can cause dramatic changes to local heat transfer because they can modify 
the production of turbulence close to the wall (discussed further in Section 2.2.1). 
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Free Surfaces: This set of technical volumes focuses on single-phase flow (see Volume 1) and 
free surfaces are therefore not considered in detail. However, the following phenomena are char- 
acteristic of buoyancy driven flows: 


« There is often a sharp temperature gradient across free surfaces. Similar to stratification 
(discussed in Section 2.1.2), if this free surface level fluctuates (e.g. due to flow unsteadiness 
caused by natural convection) components or walls crossing the free surface may suffer 
repeated temperature fluctuations that may cause thermal fatigue. This is a particular issue 
for pool-type reactors with metal coolants. 


Surface evaporation can cause significant heat transfer from heated pools (e.g. spent fuel 
pools) and this can alter the flow fields within the pool and in the gas or air above it. This heat 
transfer is therefore likely to depend on the detailed three-dimensional flow field in both the 
pool and the gas or air above it. 


« Aside from evaporation, heat transfer can vary significantly across free surfaces as a re- 
sult of the sub-surface natural convection flow field (CSNI, 1994) and sharp gradients and 
fluctuating levels may lead to deposits forming on wetted structures (surface deposition is 
considered further in Volume 2, Section 2.2.2). 


Theory 


This section provides a brief overview of the theory relevant to passive cooling systems, consid- 
ering heat transfer by natural, forced and mixed convection; buoyancy, mixing and pressure loss 
effects associated with buoyancy affected flow fields; turbulence and transition; and key aspects of 
Computational Fluid Dynamics (CFD) theory relevant to buoyancy affected flows. 


Natural, Forced and Mixed Convection 


Key Dimensionless Groups: Similarity analysis of the Navier-Stokes and energy equations 
gives three key non-dimensional groups that indicate the overall nature of convection flow fields°: 


1. Reynolds number (Re), a ratio of momentum forces to viscous forces. 
2. Grashof number (Gr), a ratio of buoyancy forces to viscous forces. 
3. Prandtl number (Pr), a ratio of momentum diffusivity to thermal diffusivity. 


For forced convection, the key dimensionless groups are Re and Pr (which can be combined to 
form the Péclet number, Pe = RePr). For natural convection, the key dimensionless groups are Gr 
and Pr (which can be combined to form the Rayleigh number, Ra = GrPr). For mixed convection, 
Gr, Pr and Re are all relevant. The Atwood number, At (a non-dimensional density difference) is 
also used to characterise the stability of flows featuring density effects like stratification. 


It is possible to introduce large uncertainties in analysis work by using information that is inap- 
propriate to the flow field under study. Identifying the significance of buoyancy and inertial effects 


Other groups may be important in different scenarios, such as the Froude number (Fr), Kutateladze number (Ku) or Eckert 
number (Ec) where there are free surfaces, two-phase or high speed flow respectively. Mass transfer equivalents also exist 
where buoyancy forces are generated by different species instead of thermal effects. For further information see Rogers 
and Mayhew (1992), Incropera et a/. (2011), Rohsenow et al. (1998), Kakag et al. (1987), Schlichting and Gersten (2017) 
and Zohuri and Fathi (2015). 
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(i.e. whether flow in a particular part of a system should be considered as natural, forced or mixed 
convection and laminar or turbulent flow) is therefore important, for example to understand: 


« What surface heat transfer correlations should be used in a system model. 


¢ What flow phenomena might be expected in a CFD model, and how the mesh should be 
constructed to capture the resulting gradients of flow variables. 


The significance of buoyancy effects can be difficult to understand without analysis, so scoping 
calculations may be needed in addition to the methods discussed in this section. 


Convective surface heat transfer is important for many passive cooling systems, and is introduced 
in some detail in Volume 2 (Section 2.1.2) alongside another key non-dimensional group: 


« Nusselt number (Nu), a ratio of convective to conductive heat transfer in a fluid at a boundary. 


Nu may be defined locally or as an area-average for a given surface (Nu), and can be used to cal- 
culate local or area-averaged Heat Transfer Coefficients (HTCs), h or h respectively (as discussed 
in Volume 2). Correlations for Nu (and Nu) are normally based on the non-dimensional numbers 
discussed above, so different correlations are likely to be appropriate for natural, forced and mixed 
convection. 


Forced and Natural Convection: Nusselt numbers for pure forced and natural convection are 
often denoted Nur and Nu, respectively. A number of correlations for Nu are available for different 
geometries in research papers and books, such as Rogers and Mayhew (1992), Incropera et al. 
(2011) and Rohsenow et al. (1998). As with any correlation, care must be taken to ensure that 
the correlation is valid for the flow and fluid that is being assessed. Particular issues for buoyancy 
affected flows include: 


« The range of non-dimensional numbers the correlation is valid for, and how similar the ge- 
ometry is to the geometry assumed in the correlation. 


* Whether the surface is heating or cooling the flow and whether a uniform heat flux or uniform 
temperature is assumed. 


¢ How similar the flow field is to the flow assumed in the correlation (normally uniform, steady, 
fully developed flow is assumed). 


It is also important that the characteristic length scales, temperatures and fluid properties used in 
the correlation are appropriate, and these aspects are discussed in Volume 2. Aspects specific to 
liquid metal and molten salt are considered in Volume 5 and Volume 6. In particular, it is noted that 
the effect of buoyancy on heat transfer for low Pr fluids (like sodium and lead) can be quite different 
to conventional fluids (Jackson, 1983, Jackson et al., 1994). 


Mixed Convection: A range of flow regimes are often possible even for a simple geometry like 
a vertical tube, and correlations are likely to depend on Re and Ra (or Gr). From this perspec- 
tive, pure forced and natural convection (complex as they may be) could be viewed as relatively 
straightforward special cases. One way of aiding clarity is to use a flow regime map, which plots 
the significance of momentum (Re) against the significance of buoyancy (Ra or Gr). Lines showing 
boundaries for application of correlations for forced, natural and mixed convection and laminar or 
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turbulent flow can then be plotted, based on where mixed convection heat transfer does not deviate 
by more than, say, 5%-10% from pure forced or free convection (see Figure 2.5). 
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Figure 2.5: Illustrative flow regime map (based on a vertical tube, after Metais and Eckert, 1964). 


Different flow regime maps are needed for different geometries (in particular, external and internal 
flows are treated differently). Situations where buoyancy effects are assisting (same direction to 
the forced flow), resisting (opposite direction to the forced flow) or transverse (perpendicular to 
the forced flow) are generally also treated separately. The impact of increasing buoyancy on forced 
convection heat transfer can also be very different for laminar or turbulent flows. Guidance on these 
aspects is provided in Rohsenow et al. (1998) and Kakag et al. (1987). A number of correlations 
are presented in Runchal (2020), Gebhart et a/. (1988) and Martynenko and Khramtsov (2005). 


Aside from graphical methods, the dimensionless group below may be helpful in determining 
whether flow should be considered as forced, natural or mixed convection: 


Gr 


Rel a ratio of buoyancy forces to momentum forces 
e 


Forced convection is then expected to dominate as Gr/Re” — 0 and natural convection as 
Gr/Re" — oo. Dimensional arguments indicate that n = 2, and this group (Gr/Re) is often re- 
ferred to as the Richardson number, Fi (particularly used in environmental fluid mechanics, Turner, 
1973). However, other values of n have been correlated with experimental results for different ge- 
ometries, inclinations, boundary conditions and fluids, and other dimensionless groups including 
Ra and Pr have also been proposed. A review of relevant literature (such as Gebhart et a/., 1988, 
Kakag et al., 1987, Jackson et a/., 1989, Rohsenow ef al., 1998) is likely to be helpful before using 
these groups to provide anything more than a rough indication. 


Mixed convection in turbulent flows can be especially complex, as differing levels of buoyancy 
force can cause marked reductions or increases in heat transfer (principally due to significant 
modification of turbulence production in the turbulent portion of the boundary layer). This situation 
in vertical pipes has received significant study (partly as a result of the application to boiler design) 
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and a review is presented in Jackson ef al. (1989). This work has been extended to passive cooling 
systems (IAEA, 2002b, Annex Part 2) and an overview is provided in Jackson (2018). 


In general (at the time of writing), detailed mixed convection guidance tends to be contained in text 
books written some time ago (like the references above) or in papers in the literature tailored to 
specific systems (for which general applicability may be unclear). Careful review of the literature 
may therefore be needed to identify guidance relevant to a specific situation under study. 


Buoyancy, Mixing and Pressure Loss 


Buoyancy provides the driving force for natural convection (Section 2.1.1). It is a ‘body force’ re- 
sulting from gravity acting on volumes of fluid with different densities. The hydrostatic pressure 
differences it generates are proportional to the product of the difference in density (Ap) and the 
height over which this difference is maintained. Accurate material properties are therefore impor- 
tant in assessing the impact of buoyancy (Section 2.3). 


The motion driven by buoyancy effects is counteracted by energy dissipation, or pressure losses 
(essentially entropy generation). These result from viscous effects associated with the mixing of 
the flow across gradients, in situations like plume boundaries, near walls, and between different 
flow features. This mixing is often dominated by complex turbulence phenomena (Section 2.2.3). 
Assessing these effects is a challenge across the whole field of thermofluid mechanics, of which 
Nuclear Thermal Hydraulics (NTH) is a part (Denton, 1993 provides internal aerodynamics con- 
text). Natural convection flow fields in particular are driven by temperature gradients that are highly 
dependent on the mixing between different areas of fluid or the distribution of heat sources and 
sinks. How this mixing is treated will depend on the scenario under study. 


Pools and Plena: For free plumes rising through a quiescent fluid (such as in a pool), correla- 
tions may be used to predict the development of plumes (for example Turner, 1973; Shabbir and 
George, 1994; Gebhart et a/., 1988; Woods, 2010). These correlations effectively predict the bal- 
ance between the driving buoyancy effects and energy dissipation associated with mixing at the 
plume boundaries, and are often based on experimental observations in idealised situations. 


However, the situation in pools and plena in a cooling system may differ significantly from the 
idealised situation used to develop these correlations, particularly if there are interactions between 
different flow phenomena or between the flow and structures. As a result, the flows in pools and 
plena are often investigated using experiments or CFD (Section 3) to gain detailed information and 
understanding of the flow field. In some situations (e.g. for a specific plena), experiments or CFD 
may be used to develop application-specific correlations. 


Loops and Channels: Open or closed loop systems are often assessed using a system analysis 
(Section 3.1), where coefficients are used to predict the stagnation pressure loss over a part of the 
system. The two most generally useful coefficients are the ‘stagnation pressure loss coefficient’ or 
‘loss coefficient’ (K) and the related ‘flow resistance’ (¢): 


w2 
APr = K5pU? =¢" where ¢= 


2A2, 
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The loss coefficient (K) is dimensionless. K is widely used because for low-speed turbulent flows 
(i.e. the majority of flows of interest to engineers) the stagnation pressure losses (and associated 
energy dissipation and entropy generation) scale with 5pU", the ‘dynamic pressure’, which is the 
kinetic energy per unit volume of the flow. Loss coefficients are generally based on experimental 
observations in idealised scenarios, and are available for flow through a wide variety of components 
(Miller, 2009; Idelcik and Ginevskil, 2007; IAEA, 2001). 


The flow resistance (¢) is not dimensionless (it has units m~*). However, ¢ is easier to work with 
in a systems analysis, because it can be summed directly to calculate an overall pressure loss 
through a number of components, using a common mass flow rate (W). This mass flow rate is 
generally of interest in a system analysis because it is conserved through continuity. As a result, 
calculating pressure losses from K involves continually using some characteristic area (A,,) for 
each component to calculate an appropriate local velocity (generally a superficial velocity, Volume 
2, Section 2.1.2), which is inconvenient and causes a risk of incorrect/inconsistent areas being 
used. 


A number of other coefficients may also be encountered for calculating pressure losses: 


¢ For lengths of pipework or ductwork, K is normally calculated from the pipe (or ‘Darcy’, or 
‘Moody’) friction coefficient (f, introduced in Volume 2, Section 3.4.2.5), using the length to 
diameter ratio of the pipe (i.e. K = fL/D). The ‘hydraulic diameter’ is useful for non-circular 
geometries (i.e. Dh = 4Acs/pcs, where pcs is the perimeter of the cross-section). 


« For nozzles and orifices, a ‘discharge coefficient’ (C,) is sometimes used, which is the ratio 
of the actual discharge to the theoretical discharge. This can normally be converted using 
K = 1/C3 if K and Cg are based on the nozzle/orifice area or K = (Ai /A2Cq)? if K is based 
on the upstream pipe area (A) and C, is based on the nozzle/orifice area (Ap). 


* For valves, a number of dimensional ‘flow coefficients’ are in use (typically C,, imperial, and 
K,, metric). A number of other coefficients may be used to provide additional corrections for 
compressible flows such as for pressure relief valves. Pressure loss data and guidance for 
using this is normally available from manufacturers, but care must be taken to ensure that 
these are used appropriately (additional detail on this topic is available in Miller, 2009). 


Correlations for K are often based on uniform inflows, perfectly mixed outflows, isothermal and 
incompressible flow and are normally only valid for specific ranges of Re and Pr (although correc- 
tions may be available to extend their use, Miller, 2009). As such, for complex flows (like natural or 
mixed convection) their applicability may be uncertain, and a correlation may not be available for 
complex geometry in part of a system. In parts of a system where these aspects are a concern, 
experiments or CFD may be used to investigate the flow in more detail and provide comparative 
data, or develop application-specific correlations (NSC, 2018). 


Example: To illustrate balancing buoyancy and pressure loss in an air chimney, consider a Reac- 
tor Vessel Auxiliary Cooling System (RVACS) driven by a hydrostatic pressure difference between 
the heat addition from a reactor (at height h,) and the top of the chimney (at height hz). The air 
in the outer annulus is assumed to be at a cold ambient temperature, while the air in the inner 
annulus is heated so that the difference in density between the air at the inlet and outlet of the 
RVACS is Ap = p1 — po (Figure 2.6). In this simplified example, the driving buoyancy is assumed 
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to be balanced by the total pressure loss for the flow through the inner and outer annuli AP;, so 
that: 


kK, 1 Ky 1 
Apg(ho — hi) = APr1 + APr2 = Ki3piU2 + KoipoUs = | —,— +5 -— |W’ 
pg(he 1) T1a+ T,2 15P1U, + K25p2U5 2A a 2A rs 
Where g is acceleration due to gravity and U is the flow velocity through each annulus (calculated 
from the flow rate and cross-sectional area), but writing the equation in terms of a common W 


shows that there is only a single flow rate variable that needs to be calculated. 


Reactor 


eiieie é| i 


Figure 2.6: Passive ventilation (natural circulation) in a simple RVACS. 


he 


Turbulence and Transition 


Turbulence and transition can have a significant impact on the flow in passive cooling systems, par- 
ticularly where surface heat transfer occurs. The nature of turbulence and the transition between 
laminar and turbulent flows is a large and complex topic that has been a subject of international re- 
search for many decades. This section therefore provides an overview of the aspects of turbulence 
and transition most relevant to natural convection. 


Turbulence: Turbulence is a flow phenomenon characterised by apparently chaotic and irregular 
variations in flow velocity and pressure. From a physical perspective, a flow will transition from 
laminar to turbulent when its kinetic energy exceeds a level where viscosity can no longer damp 
out small perturbations arising in the flow. These perturbations can then grow, causing eddies over 
a wide range of temporal and spatial scales (Figure 2.7). 


In a buoyant free plume, the initially laminar flow may quickly become turbulent, with eddying motion 
at the plume boundaries causing significant entrainment (or mixing) of external fluid into the plume, 
increasing the mass of fluid moving in the plume but reducing its temperature and hence the driving 
body force (Section 2.1.1). Turbulence enhances the mixing of momentum and thermal energy 
(across velocity and temperature gradients respectively) thereby exerting a strong influence on 
the behaviour of flows near walls, particularly in boundary layers, wall-bounded shear layers, or 
separations. This can affect surface heat transfer (Section 2.1.3), flow pressure losses and mass 
flow rates (Section 2.2.2). This is a significant topic in its own right, and a detailed discussion on 
boundary layer flows is provided by Schlichting and Gersten (2017). In a turbulent mixed convection 
flow, existing shear-driven (i.e. non-buoyant) turbulence may quickly mix out warmer fluid from the 
near-wall regions before plume flow structures can form. 
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Figure 2.7: False colour Schlieren pictures showing density variations in turbulent and laminar flows. 


At the largest scales, turbulent eddies can be of the same order as the flow geometry. At the 
smallest scales, the kinetic energy of the turbulent eddies is dissipated as heat by viscosity. As 
Re increases, a widening inertial sub-range occurs between these scales (Figure 2.8). From a 
mathematical perspective, it is the non-linear nature of the momentum terms in the Navier-Stokes 
equations which gives rise to turbulence (turbulence theory is discussed in detail in Pope, 2000). 
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Figure 2.8: The kinetic energy spectrum of turbulence. 


Turbulence is prevalent in most engineering flows, and often has a significant impact on processes 
important to NTH analysis. One of its impacts is increased mixing of advected quantities (such as 
energy, momentum or species concentration). Once a flow becomes turbulent, the turbulence (and 
its associated mixing) typically becomes the dominating feature and will often play a crucial role in 
determining many other important engineering parameters, such as: 


¢ Heat transfer from surfaces. 
* Mixing of momentum and energy. 
« Pressure drops and mass flow rates. 


18 of 91 


Technical Context 


* Forces acting on immersed bodies and surfaces. 
* Existence and development of secondary flows or stratification. 
* Solid deposition and transport. 


Because the nature of laminar and turbulent flow is very different, it is important to understand 
whether flow in a particular part of a system should be considered as laminar, turbulent or tran- 
sitional (in a similar way to identifying whether a flow should be considered as natural, forced or 
mixed convection). While general engineering flows are often turbulent, the low velocities charac- 
teristic of natural circulation make laminar or transitional flows more common than would be normal 
for forced circulation. 


The treatment of turbulence in modelling work varies with the approach employed. For a system- 
level analysis (Section 3.1), this is likely to include understanding whether the flow is likely to be 
laminar or turbulent, steady-state or time-dependent, developing or fully developed®, so that heat 
transfer correlations can be used appropriately (Section 2.1.3). For CFD analysis (Section 3.2) this 
is rather more involved, and is discussed in Section 2.2.4. 


Transition: Since key parameters like surface heat transfer and pressure loss coefficients vary 
significantly between laminar and turbulent flows (Volume 2, Incropera et a/., 2011), transition is of 
interest for passive cooling systems. There are a number of different laminar to turbulent transition 
mechanisms, which tend to occur in different situations (Schlichting and Gersten, 2017). 


« Natural transition results from the amplification of hydrodynamic instabilities in a stable lam- 
inar flow. 


¢ Bypass transition often occurs in turbomachinery, and results from disturbances caused by 
high upstream turbulence levels causing ‘spots’ of turbulence that trip boundary layers (‘by- 
passing’ natural transition). 


* Separated transition results from separation of laminar boundary layers under adverse pres- 
sure gradients causing a free-shear layer that may (or may not) reattach as a turbulent bound- 
ary layer. 


Depending on the type of flow, it may be possible to identify approximate conditions around which 
transition might be expected to occur: 


« For forced convection where buoyancy forces are negligible, Re is the key non-dimensional 
group. Transition is expected at around Re ~ 2,300 for steady and fully developed flow in 
pipes and around Re ~ 500,000 for steady flow over flat plates at zero incidence (Schlichting 
and Gersten, 2017). It is important that Re is based on appropriate fluid properties (Sec- 
tion 2.3) and an appropriate characteristic length (e.g. internal diameter for fully developed 
internal pipe flow or length from the leading edge for external flow over a plate). Real plant 
geometries and flows may differ substantially from uniform flow in pipes or over flat plates, 
and this may increase or decrease the Re at transition. Increased surface roughness or up- 
stream turbulence can substantially reduce the Re at which transition occurs. 


For natural convection, Gr and Ra are the key non-dimensional groups. Gebhart et al. (1988) 
provides specific non-dimensional groups based on Gr for predicting the start and end of 


§ Although the correlations used in system codes may well be for developed flow only. 
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transition (Gr is also considered more useful than Ra for predicting transition in fluids where 
Pr is not close to unity such as metals, Dzodzo, 2018). However, as a rough indication in 
common fluids like air and water, transition is expected at around Ra ~ 10° for a vertical 
plate (Incropera et a/., 2011). As noted above, it is important that non-dimensional groups 
are based on appropriate fluid properties and characteristic length (e.g. the height of the 
plate). Complex geometry and surface roughness may change transition behaviour signifi- 
cantly (similar to forced convection). 


For mixed convection, Re and Gr are important and the situation is much more complex, even 
for a simple geometry. Flow regime maps may be helpful and are discussed in Section 2.2.1, 
and a review is presented in Runchal (2020). 


It is noted that these are typical values, the actual values of Re and Gr at transition may vary. In 
particular, the exact onset of transition from laminar to turbulent flow is affected by factors such as 
local geometry, surface finish, upstream turbulence, buoyancy effects and vibration. Complex ge- 
ometries with non-uniform flows, rough surfaces and high upstream turbulence are likely to cause 
transition at lower Re, and vice-versa. Also, even at high bulk Gr values, regions with stagnant fluid, 
laminar circulation, transitional and fully turbulent regions may exist in complex geometries. 


In buoyancy driven flows, laminar, transitional and fully turbulent regions can often be present 
simultaneously, as a result of the relatively low driving forces. For example, in a simple differentially 
heated square cavity, significant amounts of turbulence are generated close to the vertical walls 
(where temperature gradients are large), which is convected around the cavity by the convection 
cell that develops (driven by the mean buoyancy force). However, in the centre of the cavity the 
flow may remain relatively quiescent and therefore be laminar (Hanjalié, 2002). 


Buoyancy Affected Flows and CFD 


Buoyancy affected flows present particular challenges for CFD modelling, particularly if they are 
turbulent. The heat transfer (whether by natural, forced or mixed convection), buoyancy effects, 
mixing and pressure losses discussed in the previous sections can all be coupled within a flow 
field, and their detailed impacts on this flow field can be strongly affected by turbulence. As such, 
the prediction of heat transfer, buoyancy, mixing and pressure losses is strongly linked with the 
approach taken to predict turbulence, one of the key challenges for CFD modelling work. 


The common approaches for predicting turbulence are introduced in Volume 1 (Section 4.5.3). 
The practical aspects of using CFD to model buoyant flows are discussed in Section 3.2. By con- 
trast, this section presents an overview of selected theory on Large Eddy Simulation (LES) and 
Reynolds-Averaged Navier-Stokes (RANS), the most commonly used CFD approaches, with a fo- 
cus on heat transfer and buoyant flows. 


LES 


In LES, the larger length scales are resolved and the effects of the unresolved smaller motions are 
modelled using a Sub-Grid-Scale (SGS) model (Runchal, 2020, Ciofalo, 1994). Early SGS models 
were developed with the principle that turbulence is more isotropic at the smaller scales (Pope, 
2000), meaning that these were much simpler than their RANS counterparts. However, as with 
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RANS turbulence modelling, a number of higher-order models have also been developed. Some 
of the more commonly available models include: 


Smagorinsky: First introduced by Smagorinsky (1963), this expresses the eddy-viscosity as a 
yi/3 


cell 


function of a length scale (the effective grid size, usually computed as A = )anda 
velocity scale (the filtered strain-rate). A coefficient, termed the Smagorinsky constant (or 
coefficient), is defined by assuming local equilibrium between sub-grid turbulent kinetic en- 
ergy production and dissipation. The model does not therefore properly account for non- 
equilibrium effects (such as may occur in buoyancy driven flows). The Smagorinsky constant 
is normally set by the user across the whole solution, so cannot be tailored for local flow 
regimes. Modifications are required to make the model valid for near-wall regions. Buoyancy 


extensions have been proposed by Eidson (1985). 


Dynamic Smagorinsky: One key issue with the Smagorinsky model is that the appropriate value 
for the Smagorinsky constant is different in different flow regimes (Pope, 2000). The Dynamic 
Smagorinsky model provides a means for determining an appropriate local value for the 
constant (or coefficient), as discussed in Lilly (1992) and Germano et al. (1991). 


WALE: The Wall-Adapting Local Eddy-Viscosity (WALE) model extends the Smagorinsky model 
by relating the turbulent viscosity to the local rotation rate, in addition to the local strain rate 
(Nicoud and Ducros, 1999). It is designed to return the correct near-wall behaviour for wall 
bounded flows and is thus generally recommended over the Dynamic Smagorinsky model. 


Whether or not buoyant extensions to the above SGS models are required will depend on both how 
dominant the buoyant effects are in a particular flow, and how fine the LES mesh is. With coarser 
meshes, the contribution of the unresolved scales to the solution increases and thus the choice of 
SGS model can be more important. 


For heated flows, a term representing unresolved temperature variations appears in the filtered 
energy equation and is usually referred to as the SGS heat flux. As with RANS (Section 2.2.4.4), 
models for the SGS heat flux are less mature with most approaches using a turbulent (or eddy) 
diffusivity approach along with a fixed turbulent Prandtl number. 


RANS Context 


In RANS, the mean (Reynolds-averaged) flow is resolved, and the effects of turbulence on this 
mean flow is modelled using a turbulence model. For flows where buoyancy or heat transfer play 
a dominant role, such models must account for both the turbulent transfer of momentum (using 
turbulence models to predict the Reynolds stresses) and the turbulent transfer of heat (using turbu- 
lent heat transfer models to predict the turbulent heat fluxes). For natural convection flow fields in 
particular, these are closely coupled and reliable predictions of both are generally important. While 
a large number of models are available for turbulent momentum transfer, far fewer are available for 
the turbulent transfer of heat or other scalars (Hanjalié, 2002). 


Challenges associated with buoyancy and the turbulent heat fluxes are introduced below. Sec- 
tions 2.2.4.3 and 2.2.4.4 then provide an overview of the turbulence and turbulent heat transfer 
models typically available in commercial CFD codes, considering their limitations and application 
to passive cooling flows. It is noted that turbulence modelling is a large and active research topic 
and comprehensive coverage is not possible in this document (more information is available in, for 
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example, Durbin, 2018, Runchal, 2020, Gatski and Rumsey, 2002 and Wilcox, 2006). Near-wall 
modelling is considered in Section 3.2.6. 


Buoyancy: Buoyancy driven flows are difficult for turbulence models because of the following 
challenges (which are not normally present in non-buoyant flows): 


1. Anisotropy: Buoyancy only directly affects momentum in the direction of gravity, and has a 
similarly anisotropic effect on the Reynolds stresses and turbulent heat fluxes. The most 
commonly used turbulence models, linear Eddy Viscosity Models (EVMs), cannot directly 
account for this anisotropy. 


2. Generation and suppression of turbulence: Buoyancy forces can either generate or suppress 
turbulence depending on the alignment between gravity and the density (temperature) gradi- 
ents. Even in geometrically simple systems, this can result in a multitude of regimes (laminar, 
turbulent and transitional) coexisting, which is challenging (Hanjali¢é, 2002). 


3. Coupling: The close two-way coupling buoyancy imposes between the mean velocity and 
temperature fields also applies to the turbulence fields, so the Reynolds stresses appear in 
equations and models for the turbulent heat fluxes and vice versa. This complicates mod- 
elling and can increase the numerical stiffness of the resulting models. 


4. Multiple timescales: The presence of velocity and temperature fields gives two associated 
timescales in the flow, causing concurrent effects on both fields that also influences the 
turbulence. 


Accounting for buoyancy in RANS turbulence models typically involves Reynolds-averaging the 
buoyancy force terms and including them in the derivation of the model’s transport equations. 
The main modelling effort is predicting the turbulent heat fluxes (see below). Additional effects 
arise through the influence of buoyancy on the fluctuating pressure, which can be accounted for 
in Reynolds Stress Models (RSMs) with modifications to the so-called pressure-strain term (Sec- 
tion 2.2.4.3). 


Turbulent Heat Transfer: Turbulent heat transfer arises due to the significantly increased mixing 
that occurs once a flow becomes turbulent. In a RANS approach, the impact of this on the mean 
temperature field is captured by the turbulent heat flux. This arises from Reynolds-averaging the 
energy equation, and represents the averaged product of the velocity and temperature fluctua- 
tions. It describes the effect of turbulence on the transport of mean energy within the flow, and 
therefore affects the mean temperature in a similar manner to how the Reynolds stresses affect 
mean momentum. 


RANS Turbulence Models 


This section considers the theory associated with the main turbulence models, which are illustrated 
in Figure 2.9 (selecting models is considered in Section 3.2.6). 
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Figure 2.9: Summary of the main RANS turbulence models. 


Linear EVMs, General Remarks: These use a simple linear algebraic relationship to relate the 
Reynolds stresses to the mean strains via a turbulent (or eddy) viscosity (u;)’. They are gen- 
erally categorised into zero-equation, one-equation, or two-equation models (depending on how 
many transport equations they include). The two-equation k - € and k - w family of models and one- 
equation Spalart-Allmaras model are popular linear EVMs and are discussed in following sections. 


As a result of the linear algebraic stress-strain relationship, linear EVMs cannot predict the highly 
anisotropic impact of buoyancy forces on the Reynolds stresses (or interaction between the Reynolds 
stresses and turbulent heat fluxes). These models are therefore challenged by natural or mixed 
convection flows where the direct influences of buoyancy on shear stresses are significant. They 
are also challenged by complex flows that may occur in natural convection (such as impinge- 
ment, streamline curvature, rotation, strong unsteadiness, turobulence-driven secondary flows) be- 
cause the linear relationship is unable to fully capture interactions between the mean strain field 
and anisotropic Reynolds stresses. While the aggregate effect of these influences may appear in 
the turbulent kinetic energy equation, it will not normally capture the full effects on the Reynolds 
stresses. Ad-hoc corrections to address these aspects are ultimately empirical and will likely only 
be of benefit within the intended range of applicability. 


Despite the above well-known weaknesses, linear EVMs may perform surprisingly well in buoy- 
ancy affected flows, since they are capable of responding sufficiently to changes in mean-strain 
caused by the action of buoyancy forces (rather than the direct impact of buoyancy on the turbu- 
lence itself). However, their performance may well be poor where the direct influence of buoyancy 
on turbulence is a dominant feature (Section 2.2.1), or where buoyancy (or other) forces lead to 
strongly anisotropic strain fields. 


The inclusion of buoyancy effects in linear EVMs varies, but generally involves using terms to rep- 
resent the generation or dissipation of turbulent energy due to the direct action of buoyancy. For 
models which solve a transport equation for the turbulent kinetic energy (k), a buoyant generation 
term is usually added to the equation. This is exact, but based on the predicted turbulent heat 
fluxes. Whilst this represents the direct effects of buoyancy on the turbulent kinetic energy, the ef- 
fects of buoyancy on the anisotropy of the turbulence cannot be included. For two-equation models 
which solve a length scale determining transport equation (e.g. ¢ or w), an equivalent term can 


The terms pz and 1; are both used in this context and may both be referred to as the turbulent (or eddy) viscosity. However, 
since v is the kinematic viscosity, v; is often referred to as the kinematic eddy viscosity. 


23 of 91 


Technical Context 


be added, but may not be, since the direct effect of buoyancy here is less well understood. CFD 
software may include this by default, or include an option to, and the constants used may vary. 


Linear EVMs, Zero Equation Models: These use very simple algebraic relations to predict the 
turbulent viscosity, mostly using constants input by users or simple measures such as wall distance. 
Examples include Prandtl mixing length, Cebeci-Smith and Baldwin-Lomax model. Increased com- 
puter power and availability of more advanced models mean these models are rarely used. 


Linear EVMs, Spalart-Allmaras Model: Proposed by Spalart and Allmaras (1992) for flows over 
wings, this is a one-equation model that solves a transport equation for the turbulent viscosity di- 
rectly. This model has been shown to give good results for boundary layers in external aerody- 
namics, but has not been calibrated for general industrial flows and there are few examples of 
application to buoyancy driven flows in research literature. As such, for passive cooling applica- 
tions, its use is likely to be limited to starting up more complex calculations with stability problems. 


Linear EVMs, k-€ Models: Proposed by Jones and Launder (1972), this is one of the best- 
known and most widely used linear EVM. This two-equation model solves transport equations for 
the turbulent kinetic energy, k, and the dissipation rate of that energy, ¢. The original model did not 
account for viscous effects found in the near-wall region (the flow was assumed to be fully turbulent 
everywhere) but a number of low-Re variants have been proposed which do include terms to allow 
the model to be applied in such regions® (Chien, 1982; Launder and Sharma, 1974; Yang and Shih, 
1993). The original model did not include direct buoyancy effects, but this can be included as noted 
in the general remarks above. Many model variants exist, but two are most common’: 


1. Realizable k-« model includes an alternative formulation for the turbulent viscosity (using 
a variable c,,) and a modified transport equation for ¢ (Shih et a/., 1995). While realisability 
appears desirable, turbulent stresses are not really generated in the way an EVM supposes, 
limiting its efficacy (Hanjali¢ and Launder, 2011). 


2. RNG k-e model applies a statistical technique called Re-Normalisation Group (RNG) theory 
to the Navier-Stokes equations (Yakhot et a/., 1992). The resulting model has the same form 
as standard k - € but has coefficients arising from the procedure and an additional term in the 
€ equation designed to improve sensitivity in flows with strong straining. It is argued that this 
should improve performance over the standard k - € for a wide range of flows. 


One of the weaknesses of the k- € family of models is a tendency for large turbulent length scales 
to arise in strong adverse pressure gradients, linked to modelling of diffusion processes within the € 
equation'®. However, this is addressed by the Yap correction (Yap, 1987; Craft et al., 2000), which 
can significantly improve velocity profile predictions in buoyancy-driven turbulent flows in vertical 
pipes (Ince and Launder, 1989) and should be considered if available. 


Often references within literature refer to a ‘standard’ k -¢ model; this typically refers to the high-Re version of Jones and 
Launder (1972), but many studies incorrectly associate the term ‘standard’ with later variants. This is important because 
many deficiencies of the ‘standard’ model have been addressed in later variants. 

As published, these do not include wall-damping terms and thus must be used with wall functions, but some CFD developers 
may include modifications to enable such models to be used with a low-Re near-wall approach, so software documentation 
should be checked. 

10 For impinging jets this was shown to result in overly excessive heat transfer rates (Craft et al., 1993) and termed the 
turbulent/round jet anomaly. 


o 
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Linear EVMs, k-w Models: Proposed by Wilcox (1988), this is another widely used family of 
linear EVMs. This two-equation model solves an equation for k and w (instead of ¢) and can solve 
right up to the wall without further modifications (i.e. without introducing viscous damping terms, as 
must be done for the ¢ equation). As for k - €, direct buoyancy effects can be included as noted in the 
general remarks above. The original model displayed strong sensitivity to free-stream conditions 
(Menter, 1994), but this has been addressed in later variants: 


1. Baseline (BSL), Menter (1994) effectively blends a k-w formulation in the near-wall region 
with a k -€ model in the far field, as this does not demonstrate such free stream sensitivity. 


2. Shear Stress Transport (SST), Menter (1994) further developed the BSL model by introduc- 
ing a limiter into the turbulent viscosity formulation, effectively causing the model to respond 
more like a shear-stress transport equation. 


3. Generalized k -w (GEKO) model, Menter et a/. (2020) replaces the relatively inflexible cali- 
brated coefficients used by most models with a set of user-definable parameters. This allows 
the model to be adjusted in specific areas without invalidating the base calibration, and at- 
tempts to provide greater flexibility within one model to cover a wider range of applications. 


At time of writing, k-w models have been less well used for buoyancy driven flows than k-«e 
models, and the implementation of buoyancy effects in commercial CFD software has been incon- 
sistent. For some ascending mixed convection flows in particular (Keshmiri et a/., 2012), the k-w 
SST model has struggled to capture the impairment in heat transfer experienced as the influence 
of buoyancy was increased. Although interesting, GEKO is relatively new and there are therefore 
few published examples of its application, especially relating to buoyancy affected flows. 


Other Linear EVMs: A range of other linear EVMs have been developed, of which a number are 
based on the elliptic relaxation concept introduced by Durbin (1991) and one of the more popular is 
the v?-f model. In addition to k and ¢, this solves transport equations"! for v2 and f, which come 
from the solution of an elliptic equation. This has been extended to account for buoyancy effects 
(Kenjeres et a/., 2005) and has demonstrated improvements over other linear EVMs in a number of 
flows involving turbulent heat transfer, including mixed convection vertical channel flow (Keshmiri 
et al., 2012), impinging jets (Behnia et al., 1999) and ribbed walls (Manceau et al., 2000). 


Non-Linear EVMs: Non-linear EVMs have also been developed to link the mean strains and 
Reynolds stresses more reliably than linear EVMs, whilst retaining the same convenient modelling 
framework. Most schemes extend the linear Boussinesq stress-strain relation to include non-linear 
(i.e. quadratic or cubic) combinations of mean strain terms. The main advantage of this is that the 
impact of flow curvature, rotation, or buoyancy on Reynolds stresses can be accounted for at least 
qualitatively, with little impact on computational cost. 


Improvements over linear EVMs have been demonstrated in flows with curvature, separation, tran- 
sition, impinging jets and anisotropy-driven secondary flows (Craft et a/., 2000, Jaramillo ef al., 
2008). Anisotropy-driven secondary flows can alter heat transfer and increase mixing in fuel bundle 
channels (Ninokata et al., 2009; EPRI, 2014). Non-linear EVMs were therefore selected by some 
participants in a CFD benchmark of the NESTOR and OMEGA 5x5 fuel bundle experiment (EPRI, 


Nyisa surrogate scalar which is the wall-normal stress component in a turbulent boundary layer, but in an elliptic relaxation 
EVM it provides an additional velocity scale (Hanjali¢é and Launder, 2011). 
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2014), and were recommended as robust and computational efficient (Kang and Hassan, 2016). 
As such, while non-linear EVMs are not extensively used for buoyancy driven flows in industry, they 
may provide improvements in future. 


Reynolds Stress Models: Unlike EVMs, RSMs!* solve transport equations for each of the 
Reynolds stresses directly, without making use of the Boussinesq eddy-viscosity hypothesis. This 
enables anisotropic turbulence effects (which may be significant in buoyant or near-wall flows, Om- 
ranian et al., 2014) to be modelled directly. This enables better prediction of the different rates at 
which individual stresses are generated, convected, diffused and dissipated (Hanjali¢é and Launder, 
2011). The performance of more advanced turbulence heat flux models are also likely to improve, 
since many of these models are reliant on accurate values of the Reynolds stresses being available 
(Section 2.2.4.4). RSMs therefore have a much wider range of applicability than EVMs. 


The direct effects of buoyancy are included via a generation term, which is exact, but (as with 
EVMs) requires accurate values of the turbulent heat flux. Most modelling effort is focused on the 
pressure-strain term, which represents a physical process absent from models using a transport 
equation for k, and CFD tools may distinguish RSMs by how this term is treated. Notable models 
include: 


1. The Launder-Reece-Rodi (LRR) (Launder et a/., 1975) and Speziale-Sarkar-Gatski (SSG) 
(Speziale et al., 1991) RSM models use linear and quadratic relations for the pressure-strain 
correlation respectively. While initially developed as high-Re models (like k-¢), low-Re ex- 
tensions are available (Shima, 1998; Hanjalié and Jakirli€é, 1998). 


2. Elliptic-Blending RSM (EBRSM) (Manceau and Hanjali¢é, 2002) was developed based on 
an elliptic relaxation concept (Durbin, 1991) and uses a different approach to modelling the 
pressure-strain contribution. This has been successful in a number of complex flows (Billard 
et al., 2012), including those affected strongly by buoyancy (Dehoux et al., 2017). 


3. Two-Component Limit (TCL) (Craft, 1998) is a low-Re form of LRR developed to satisfy the 
so-called TCL, a highly anisotropic state where one of the normal stress components be- 
comes very small compared to the other two.'? This has demonstrated superiority in a num- 
ber of buoyancy related flows (Craft et a/., 1996; Omranian et al., 2014). 


In addition to solving equations for each Reynolds stress, RSMs also need to solve a length scale 
determining equation. Typically either ¢ or w is used, inheriting the deficiencies present in those 
equations. From a physical perspective, RSMs are particularly appropriate for modelling flows with 
strong buoyancy influences, especially in unfamiliar situations. A detailed assessment of RSMs 
can be found in Hanjalié and Launder (2011). 


12 Also known as ‘stress transport models’ or ‘second-moment closures’. 
13 While appearing academic, flow in the very near-wall region may well approach this limit, as may a strongly stable buoyancy 
affected flow. 
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RANS Turbulent Heat Transfer Models 


Modelling turbulent heat transport has received much less attention than turbulent momentum 
transport, and while a number of approaches exist, options in most CFD tools are limited. As for 
turbulence models, there is a hierarchy of models (Figure 2.10) ranging from simple turbulent (or 
eddy) diffusivity approaches (relating turbulent heat fluxes to mean temperature gradients, analo- 
gous to the turbulent viscosity concept) to full differential transport approaches (solving separate 
equations for each component of the turbulent heat flux vector). The choice of turbulent heat trans- 
fer model is an area of active research (Section 4), and is related to the choice of turbulence model 
(Section 3.2.6). 


Linear EVM =enn 
(k-e€, k-w) 


Turbulent Momentum Transport 
yodsues| }eayH JUsiNquNy 
JO uolNjosey Bulseasou| 


Increasing Resolution of 


Figure 2.10: Summary of the main RANS turbulent heat transfer models. 


Simple Gradient Diffusion Hypothesis (SGDH): This relates the turbulent heat flux compo- 
nents to the mean temperature gradients using an isotropic turbulent diffusivity (normally assumed 
proportional to the turbulent viscosity via the turbulent Prandtl number, Pr;)'*. Pr; may be a con- 
stant by default, but CFD tools often allow this to be varied by the user (Pr; is not constant in many 
buoyancy affected flows). This approach assumes that the turbulent heat fluxes are not directly 
influenced by buoyancy, and depend directly on the corresponding mean temperature gradient (a 
constant temperature in one direction implies no turbulent heat transport in that direction), which 
may well not be the case. This model therefore has limited validity, and while it may work well 
in some forced convection flows (e.g. ribbed ducts/passages, Manceau et a/., 2000) it can cause 
errors in flows dominated by buoyancy (Hanjali¢é, 2002). 


Generalized Gradient Diffusion Hypothesis (GGDH): Introduced by Daly and Harlow (1970), 
this generalises the SGDH by introducing an anisotropic turbulent diffusivity (which depends on all 
components of the Reynolds stress tensor). This model therefore depends on reliable predictions 
of the Reynolds stresses. It considers all of the mean temperature gradients (rather than just the 
one aligned with the turbulent heat flux), but if these mean gradients are all zero, it will still re- 
turn a zero value. It has little computational overhead compared to SGDH, and has demonstrated 
improvements in natural convection with strong stratification (Ince and Launder, 1989). 


14 An application of the Reynolds analogy, which assumes similarity between the way turbulence affects the transport of 
momentum and heat. Similar models may be available for mass transfer using the turbulent Schmidt number, Sc; 
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Algebraic Heat Flux Model (AHFM): AHFM is based on an algebraic simplification of the full 
transport equations for the turbulent heat fluxes, and retains all of the major generation terms 
(mean temperature gradients, mean velocity gradients and direct buoyancy effects). Since the 
SGDH and GGDH only include contributions due to mean temperature gradients, AHFM can poten- 
tially account for a wider range of buoyant interactions. The model requires values of temperature 
variance, which can be calculated algebraically (similar to GGDH) or, more commonly, by solving a 
transport equation for it. AHFM has been successfully applied to a range of buoyancy driven flows 
in enclosures (Hanjali¢, 2002; Hanjalié et a/., 1996). 


Differential Flux Model (DFM): Solves transport equations for each of the turbulent heat fluxes 
(and normally the temperature variance). This better accounts for the various ways turbulent heat 
fluxes are generated, transported and dissipated and addresses a number of deficiencies with 
more primitive models (SGDH, GGDH and AHFM). However, DFM adds significant computational 
cost, and is not normally available in commercial CFD software at time of writing (Dehoux et al., 
2017). 


Fluid Material Properties 


Natural convection flow fields are driven by density variations in response to temperature changes, 
so the material properties of the fluid (density, viscosity, specific heat capacity, thermal conductivity) 
are particularly important. In general, the properties of fluids vary due to a range of factors that 
should be considered, including: 


Nature of Fluid: The properties of fluids may be altered by a range of factors, including contami- 
nation, addition of different species, specific mixture/solution composition and ageing under 
radioactivity. The key factors to consider, and their impact, are likely to be case specific. 


Pressure: There are often two aspects to pressure variations: firstly, the overall pressurisation of 
a system and secondly, pressure changes occurring in the flow. The impact of both of these 
on material properties should be considered. For liquids, variations in material properties 
with flow pressure changes are often small, so that the impact of flow pressure changes 
on material properties may be neglected (from a density perspective this is equivalent to 
treating the flow as incompressible). In this case, the properties can be evaluated at the 
appropriate system pressure. For gases, often only the density is significantly affected by flow 
pressure changes, and these compressibility effects usually only have a significant impact 
where Ma & 0.3 (which is unusual for natural convection). 


Temperature: Temperature changes often cause larger variations in properties than pressure 
changes. For natural convection in particular, it is precisely the variation in density with tem- 
perature that provides the buoyancy forces that drive the flow (Section 2.2.2). Other proper- 
ties such as viscosity, specific heat capacity and thermal conductivity may also exhibit signif- 
icant variation with temperature. The variations in properties with temperature are therefore 
a key aspect to consider when performing analysis of buoyancy affected flow fields. 


Critical point and phase changes: Near the critical point or phase changes, very large variations 
in properties may be encountered with small changes in pressure or temperature, and these 
can have a large impact on flow fields and heat transfer. 
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The rest of this section identifies some sources of material property data and considers how these 
might be used. Including properties in modelling work may be challenging where there are strong 
gradients in properties, such as near the critical point or phase changes as noted above. These 
challenges may include accuracy of the data itself, the accuracy of the calculation of properties 
within the model, and the stability of the model. 


The properties of solids and the impact of surface modifications on heat transfer are considered in 
Volume 2 (Section 2). 


Sources of Property Data 


Possible sources of fluid material properties data for a number of common fluids are presented 
below. Whatever source is chosen, this is a key input to analysis work, so it is necessary to record 
the approach taken and consider the accuracy of the data. Uncertainty in property data can have a 
significant impact on the analysis of natural convection flows; the management of input uncertainty 
is considered in Volume 4 (Confidence and Uncertainty). 


Water and Steam: The International Association for the Properties of Water and Steam (IAPWS) 
presents formulations for properties of light water, heavy water and sea water at the time of writing, 
and new formulations may be released periodically'>. Two of these formulations for the thermody- 
namic properties of light water (IAPWS-95, Wagner and PruB, 2002, and the less computationally 
intensive IAPWS-IF97, Wagner et a/., 2000) are considered accurate by IAEA (2006). Properties 
data based on current formulations, including saturation curves, may be available from NIST (via 
the Chemistry WebBook or REFPROP tool'®) or steam tables like Haywood (1990). Saturation 
curves for light and heavy water (and more general properties for light water) are also presented in 
IAEA (2009b). It is inappropriate to consider steam as a perfect or ideal gas. While phase change 
is not considered in this volume, it is noted that condensation can be greatly affected by the pres- 
ence of non-condensible gases (further detail available in Rohsenow et al/., 1998 and Collier and 
Thome, 1996). 


Carbon Dioxide, Nitrogen and Helium: If pressures are not high and temperatures are not 
low (or if properties are not of great concern, like basic studies, model debugging etc.) it may be 
appropriate to consider these gases as perfect or ideal'’. Real gases are more likely to behave like 
ideal gases than perfect gases; further discussion on the use of perfect or ideal gas assumptions is 
presented in Rogers and Mayhew (1992). It may well be appropriate to check perfect or ideal gas 
assumptions against full properties, or full properties data can be used exclusively. Full properties 
data for these gases may be available from NIST Chemistry Webbook (for pure gases) or NIST 
REFPROP (for pure or mixed gases). Data for specific conditions may be available in textbooks, 
such as Incropera et a/. (2011). Data for helium at low pressure is included in IAEA (2009b). These 
gases are often used in closed or sealed systems where moisture content is controlled at a low 
value, but if this is not the case, the impact of moisture on properties can be significant and should 
be considered. 


15 www.iapws.org 
16 webbook.nist.gov/chemistry and www.nist.gov/srd/refprop respectively 
17 Perfect gases obey the ideal gas law and have a constant specific heat capacity; ideal gases only obey the ideal gas law. 
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Air: Air is a mixture of gases and should generally be considered as such. Like carbon dioxide, 
nitrogen and helium, it may be appropriate to either use perfect or ideal gas assumptions, or use 
full properties data (see above). Published sources include Lemmon et al. (2000), Lemmon and 
Jacobsen (2004) and IAEA (2009b). Unlike carbon dioxide, nitrogen and helium, air is often used 
in open or unsealed systems, where moisture content is likely to vary. The presence of moisture 
can significantly affect properties and hence the performance of a cooling system (e.g. the impact 
on density is considered in Picard et al/., 2008). Moisture content is also an important aspect of the 
design of HVAC systems, which can have a role in preservation of mechanical equipment as well 
as habitability for operators (further detail in Haines and Myers, 2010 and ASHRAE, 2017). 


Liquid Metals and Molten Salts: The use of material properties for liquid metals and molten 
salts requires particular considerations, which are discussed in Volume 5, Volume 6 and (IAEA, 
2002a, 2013b). In general, obtaining reliable information about the properties of more advanced 
but less established coolants may be challenging, as a range of formulations may be available and 
the accuracy of the data may be difficult to ascertain. 


Implementing Property Data 


The exact approach taken in a given project will be a judgement specific to the analysis being 
performed and the programs being used. Sources of data are discussed in Section 2.3.1. Common 
ways of including this data in analysis are illustrated in Figure 2.11 and discussed below. 


(a) Interpolate Data (b) Fit Measured Data c) Equation of State 
(d) Boussinesq Appoximation e) Constant Values 


Figure 2.11: Common ways of including material properties data. 


Interpolate directly from measured data: This typically involves reading property data into an 
analysis program and using linear interpolation to provide property values. This gives full 
control over the values used, but may cause computational overheads and require a large 
amount of data to be read into the analysis program that has to be checked. Additionally, the 
user must be confident that they will be alerted when values outside the range of the data 
are requested, as this may cause errors. 
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Fit measured data: This typically involves reading property data into a computer program, fitting 
curves or surfaces to the data, and then entering the constants of these fits into an analysis 
program. The exact approach is likely to depend on what the analysis program will easily 
accept (polynomial, exponential fit etc). This may be especially useful if properties are con- 
sidered only a function of temperature (i.e. the flow is considered incompressible), gives full 
control over the values used with minimal computational overheads and saves reading a lot 
of data into an already complex analysis program. However, care must be taken to ensure 
that the fit closely represents the data and the user must be confident that they will be alerted 
when values outside the range of the fit are requested, as this may cause errors. 


Equation of state: An equation of state uses an algebraic expression to provide properties for 
given conditions (pressure and temperature). The complexity of equations of state varies 
from simple (such as the ideal gas law, e = P/RT) to extremely complex functions. It is 
important to understand the assumptions underlying these equations and confirm they are 
relevant to the assessment before using them. It is also important to implement them carefully 
and check their output, as large errors could result from inaccurate use. 


Boussinesq approximation: This approximation'® is sometimes used in analysis of buoyancy- 
driven flows. This allows the density to be considered constant in the governing equations, 
except when multiplied by the acceleration due to gravity in the buoyancy force terms (i.e. 
pg). The density variations caused by thermal expansion are included by assuming a linear 
dependence with temperature (e.g. Ap/p ~ —BAT). This simplifies the governing equations 
and may make including density changes easier, particularly for ideal gases where B = 1/T 
(although in general, 6 and other properties such as viscosity may be set as constant or 
varying when the Boussinesq approximation is being used). 

The Boussinesq approximation is only valid for small density variations (|BAT| < 1, Gebhart 
et al., 1988) and is therefore a potential source of uncertainty in analysis work. For flows 
in loops, using this approximation rather then IAPWS-IF97 (Section 2.3.1) has been seen 
to cause significant changes in predicted system behaviour (Krishnani and Basu, 2016). 
If this approximation is used, care should be taken to ensure that it is appropriate for all 
fluid properties across the range of local conditions (pressure, temperature etc) that may be 
encountered in the flow field. This might be achieved by comparing with measured data or an 
equation of state, using sensitivity studies, or analytical methods (Gray and Giorgini, 1976). 


Constant values: Using fixed values for fluid properties may be useful in some situations, but 
must be used carefully and is of little use for natural or mixed convection flows, where by 
definition density changes have a significant impact on the flow. 


In summary, using temperature dependent properties appropriate to the average pressure of the 
flow in the system of interest (or full pressure and temperature dependent properties) should be 
considered for natural or mixed convection assessments. Using fits to measured data may be an 
appropriate starting point, as this may give a good balance between accuracy and complexity. 


It is important that properties data is provided for the range of conditions that occur in the analysis 
to avoid unphysical model behaviour. Checks may (or may not) be included in modelling software to 
generate warnings where conditions are outside of the range for which properties are available. It 


18 This is not to be confused with the Boussinesq eddy-viscosity hypothesis, see Section 2.2.4. 
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may be necessary to use more than one dataset to obtain properties over the range of conditions, in 
which case the properties should be carefully combined in a consistent manner. Whatever is done, 
it is important to record the approach taken, as properties are key inputs to any flow analysis. 


Modelling Challenges 


This section introduces the typical challenges that may be encountered in modelling buoyancy 
affected flows. Addressing these challenges is considered in Section 3. 


Local Flow Unsteadiness and Complexity: Like the plume rising from a cup of tea, buoyancy 
driven flows are often unsteady and complex/three-dimensional, even where the boundary condi- 
tions are essentially constant and two-dimensional. Inside these flows, the dynamic field (the local 
velocities, including turbulence if present) is tightly coupled to the thermal field (the local tempera- 
tures) as well as the related shear forces, mixing, changing fluid properties and buoyancy forces. 
This places particular demands on modelling work that seeks to predict detailed flow fields (such as 
CFD). Similarly for experimental work, this may impact aspects such as choices of measurement 
approaches, data logging equipment and probe placement. 


System Flow Unsteadiness and Complexity: Just as buoyancy affected flow fields can be un- 
steady and complex, so can the behaviour of the system as a whole. In particular, systems using 
natural circulation may exhibit complex oscillatory behaviour with significant sensitivity to initial 
states and evolving boundary conditions during postulated accident scenarios. Small changes 
in one area may spread to the whole system, and differences between thermal and dynamic 
timescales may result in oscillatory and unintuitive flow behaviour. Even simple systems can exhibit 
complex flow behaviour as the overall system flow develops during start-up, particularly if pressure 
losses (Section 2.2.2), which can provide a level of damping, are low. The impact of instability on 
start-up times may be a challenge if these times are safety significant (IAEA, 2005 and Krishnani 
and Basu, 2016 consider stability of passive cooling systems). An example for a simple natural 
circulation loop is shown in Figure 2.12, where complex eddying local flows couple with the wider 
system flow field to cause the bulk flow to reverse repeatedly (Wilson, 2021). 


Pressure Losses and Mixing: As noted above and in Section 2.2.2, in buoyancy affected flows, 
pressure losses and mixing act to resist the movement of fluid via shear forces, and to reduce 
temperature gradients. Both of these effects can have a large impact on flows both locally and 
over the whole system. However, the unsteady and complex flow fields that can exist in buoyancy 
affected flows can challenge the assumptions used in the correlations in system models and the 
turbulence models often used in CFD (Section 2.2.3) which are often important in predicting these 
effects. 


Characteristic Timescales: As noted in Volume 3 (Section 2.5), there is generally a large dis- 
parity between the timescales associated with detailed flow effects (e.g. shedding or plume devel- 
opment, with timescales typically less than a second), system effects (e.g. through-flow duration 
or flow reversal), solid conduction (e.g. heating of large masses of metal, with timescales typically 
in hours) and the whole plant (e.g. transients lasting days, weeks, or even months or years for 
accident scenarios). If the flows during a long plant transient are to be predicted using a CFD 
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model, it is likely to be expensive to predict the whole plant transient using the smallest timescale 
characteristic of the flow field. 


Confidence: _ It is important to develop a level of confidence in the predictions of modelling work 
that is proportionate to its significance in the economic or safety case for the plant, in line with 
a graded approach (Volume 1, Section 2.2). Since buoyancy driven flows may be challenging to 
model, developing this confidence can be a particular challenge for these flows. This topic is con- 
sidered further in Volume 4. 
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Figure 2.12: Flow unsteadiness and system level instability in a simple natural circulation 
loop with a horizontal heater and horizontal cooler (Wilson, 2021). The mass flow rate is 
m in this figure. 
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This section presents an overview of methods used to predict and understand the performance 
of passive cooling systems, to assist engineers performing analysis. Typical approaches to NTH 
analysis are considered, before assessments using system and subchannel analysis, CFD analysis 
and experimental methods are discussed in more detail. 


The overall approach to NTH analysis is introduced and discussed in Volume 1 (Section 4.1). In 
general: 


* System codes are often used to predict the performance of whole systems or groups of 
sub-systems, while CFD is often used to understand the detailed flow field in specific plant 
areas. 


* Application-specific codes (like subchannel or containment codes) are often used where 
some flow field modelling is needed, but a full CFD approach is too cumbersome or compu- 
tationally expensive. 


* Experiments are often used to understand the overall behaviour using scaled test facilities, 
investigate detailed flow effects and provide comparative data to validate analysis work. 


As for any NTH analysis, the amount of technical justification required is likely to depend on the 
safety or economic significance of the system as part of a graded approach (Volume 1). It is also 
common, and often necessary for passive cooling systems, to employ different methods together, 
for example using a system code to predict overall performance and a subchannel or CFD code to 
predict flow fields in specific areas, i.e. multiscale analysis (Volume 2, Section 3.1). 


Pools and Plena: Buoyancy affected flows in pools and plena are often complex and three- 
dimensional because the flow is less constrained by walls. Relevant flow phenomena include 
plumes, thermal stratification, thermal striping, jet interaction and free surfaces (Section 2.1). In 
general, this tends to imply more of a role for methods that resolve spatial and temporal variations, 
particularly CFD and experimental methods, rather than system codes. A notable exception is the 
system code SAM (Volume 1, Section 4.4.1), which is being developed to model liquid metal pools. 


Loops and Channels: System and subchannel codes have been developed and validated for 
primary circuit flows under a wide range of fault scenarios, particularly for forced circulation con- 
ditions. However, the complex flow fields associated with buoyancy driven flows and natural circu- 
lation mean that three-dimensional effects become more significant and reduces the suitability of 
these tools. In particular, complex geometry or flow phenomena, such as flow reversal, cold water 
injection, thermosyphons (dead legs) and cold traps (Section 2.1.2) need to be carefully considered 
and may require more detailed CFD analysis to build confidence in the overall modelling approach. 
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System and Subchannel Analysis 


System and subchannel codes are introduced in Volume 1 (Section 4.4) and their capability and 
limitations for Conjugate Heat Transfer (CHT) modelling is discussed in Volume 2 (Section 3.2). 


System codes in particular are key design tools, allowing the performance of a system (or groups 
of sub-systems and modules) to be rapidly assessed in a wide range of scenarios. However, 
buoyancy-driven flows can present particular challenges for analysis using these codes. Subchan- 
nel codes or CFD are more commonly used for in-depth analysis of subchannels, where system 
codes may be used to provide boundary conditions, such as key core inlet parameters. 


Some challenges associated with modelling natural circulation scenarios are common to both sys- 
tem and subchannel codes. In general, modelling of natural circulation using these codes depends 
on a number of constitutive equations including those for friction factor, form loss and heat trans- 
fer. Most of these correlations are relatively well-established for fully-developed flow in a circular 
pipe. However, in natural circulation, the flow may well be complex and uneven or include counter- 
current flow, even within simple circular pipes. Discrepancies may also occur if one (or more) of the 
following conditions is involved: 


« Non-circular, complex geometries (e.g. tube bundles); 
* Complex flow paths (e.g. flow obstructions, orifices); 

* Developing flow; 

* Multi-dimensional effects; 

* Transition between different flow regimes. 


The geometry, conditions and flow phenomena within the system being modelled should be re- 
viewed against the Verification and Validation (V&V) database and user documentation for the sys- 
tem/subchannel code. If issues are identified regarding the suitability of the code for the scenario 
being modelled, then a number of options are available. 


« Experimental (or CFD) data could be generated for the expected natural circulation flow 
conditions to improve the existing correlations in the system/subchannel code, although this 
may require access to the source code. For example, CFD may be used to generate an array 
of coefficients to account for pressure loss due to flow obstructions (Avramova et al., 2016). 
This will require users to carry out additional V&V activities to ensure the code accuracy. 


* If sources of errors are identified, users may ‘tune’ user-provided parameters in the input 
file to remedy discrepancies. For example, users may ‘adjust’ parameters such as junction 
pressure loss coefficients, hydraulic diameter and heated diameter. 


« Users may also couple the system/subchannel code to CFD models in order to properly re- 
solve a portion of the system where complex flow phenomena and/or geometry are expected 
(Volume 2, Section 3.1). 
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System Analysis 


System codes are generally based on a two-fluid model (solving mass, momentum and energy 
conservation equations for gas and liquid phases) and use flow regime maps for different flows 
(such as vertical, horizontal and mixing flows, see Section 2.2.1). The flow regime map is used to 
select appropriate correlations (or closure relations or constitutive models) for interfacial heat trans- 
fer, interfacial drag, wall drag, wall heat transfer and wall friction based on the local flow conditions 
(Section 3.1.1.2). 


System codes have been primarily developed for forced circulation systems, which are dominated 
by pressure drop and inertial effects, rather than natural or mixed circulation, where buoyancy (den- 
sity differences) has a significant impact on the flow. While the capability of system codes to predict 
single-phase and two-phase natural circulation has been studied, these studies are generally not 
as extensive as for forced circulation. The approach, capability and limitations of system codes to 
predict natural circulation is discussed in the following sections. 


Nodalisation and Discretisation 


Nodalisation is the representation of plant geometry within a system code. The user is required 
to construct this nodalisation, which describes the network of plant components. Each component 
may then be further discretised using hydrodynamic cells. There is considerable scope for variabil- 
ity in nodalisation (as for meshing in CFD) and this is often a source of ‘user effects’. 


Nodalisation within system codes usually relies on the use of building blocks. These building blocks 
include pipe, junction, branch and single volume and each building block is equipped with com- 
ponent models. For example, pipe components are a one-dimensional system of volumes and 
junctions and its geometry is described using flow area, length and hydraulic diameter. 


Traditionally, the nodalisation process involves a compromise between accuracy and computational 
cost, however computational cost is less of a concern for modern studies. It may seem desirable 
to improve the accuracy of results by refining the discretisation to provide increased spatial reso- 
lution, but for most system codes this does not necessarily guarantee improved solution accuracy 
(Petruzzi and D’Auria, 2008), because: 


* Alarge number of empirical constitutive models have been developed assuming a fixed/coarse 
nodalisation. For example, in RELAP5, spatial derivative terms in the virtual mass force were 
neglected as these were inaccurately approximated using relatively coarse nodalisation. 


* System codes are often based on the first-order upwind scheme, and so reducing the spatial 
resolution below a certain threshold may lead to some unphysical instabilities. 


¢ The two-fluid model can become ill-posed, which can lead to numerical instabilities (Dinh 
et al., 2003). RELAP-7 claims to eliminate this issue by using a seven equation two-phase 
flow model. 


In addition, the complexity of a component or system geometry is often simplified so that it can 
be represented in a system code. This simplification may inaccurately predict the flow behaviour, 
particularly in a multi-dimensional flow and a scaled-down facility where system specific effects 
may become exaggerated. As there is no optimum approach to construct a nodalisation, it is rec- 
ommended that nodalisation and convergence studies are carried out. 
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Pools and Plena: In pools, multi-dimensional phenomena such as thermal stratification and mix- 
ing are expected. In order to capture these conditions, a pool may be modelled using a multi- 
dimensional component or parallel (vertical) channels connected by cross-flow junctions (Kumar, 
2017). The multi-dimensional components in system codes are based on an orthogonal, three- 
dimensional grid. In order to represent the actual amount of fluid within a cell or include the local 
geometry effect, a porosity factor and junction factor are used. This can be another potential source 
of user effects. It should be noted however that the code manuals claim that TRAC-based codes 
and RELAP cannot model recirculation flows in a large open region (INL, 2018, Roth and Aydogan, 
2014b). 


In general, plena are modelled using a branch component or 1D volumes connected by cross- 
flow junctions. However, under natural circulation flow conditions, multi-dimensional and/or non- 
uniform flow conditions may be observed in the plena (Salah, 2013). These effects may be captured 
by employing a multi-dimensional component in the system codes noting that multi-dimensional 
components are intended to capture large-scale effects only. 


If the multi-dimensional flow behaviour in pools and plena cannot be captured adequately by 
system codes, CFD is recommended and typically system codes are used to define appropriate 
boundary conditions for the CFD analysis. 


Loops and Channels: Modelling loops and channels in system codes is limited to employing a 
pipe/channel component. It should be noted that the prediction may be sensitive to the choice of 
nodalisation, and a convergence study is necessary to determine an optimal nodalisation for loops 
and channels (Hou ef a/., 2017). In general, it is recommended that the node length be set up such 
that all volumes have similar material Courant limits. As a starting point, a node length may be set 
such that the ratio of the node length to diameter is > 1. 


In the case of multiple, parallel channels, it is recommended that the distribution of node lengths 
is the same or similar between the channels, particularly in the case of vertical, parallel channels. 
This is done to prevent numerical oscillations between parallel channels (INL, 2018). In general, 
these parallel channels should be connected via a branch component at both the top and bottom. 
It should be noted that there is a limited amount of validation data for natural circulation in parallel 
channels (Saikia et a/., 2019). It is therefore important that users carry out convergence studies for 
discretisation to ensure an optimal node length. 


Systems operating under natural circulation may be less stable compared to forced circulation 
systems (Section 2.4) and can exhibit flow instabilities and reversal. One of the more prominent 
experiments by Welander (1967) demonstrated that flow instabilities could occur in a single-phase 
natural circulation loop. Some examples of the studies that have been undertaken to test whether 
system codes can predict these flow instabilities include: 


« Ferreri and Ambrosini (1999) showed that RELAP can predict flow instabilities however, the 
nodalisation and time step choices can lead to significant dampening of the instabilities seen. 
For the same time step, a different nodalisation scheme can lead to significant differences 
in the flow rate predicted in the pipework due to numerical errors associated with the solver 
used in the system code. It also contains advice on how to determine whether the model will 
predict the flow instability. 
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* The MTT-1 natural circulation loop (rectangular loop with heating and cooling sections at 
the University of Genoa) was simulated in both RELAP and CATHARE. Misale et al. (1999) 
found that the CATHARE model could predict the flow and temperature behaviour for heating 
powers below the stability threshold to a reasonable accuracy, but did not predict the flow 
instabilities seen in the experiment when the heating power exceeded the minimum thresh- 
old, while the RELAP model predicted flow instabilities at all heating powers, even when the 
heating power was below the threshold. This demonstrates that different codes can show 
different behaviours for similar scenarios as the instability is related to the solver used by 
the system code, although as noted previously, the nodalisation and time step can have a 
significant impact on the solution. 


Vijayan et al. (1995) used ATHLET to replicate natural circulation experiments with varying 
heating power for three different pipework diameters. The analysis involved modelling the 
experiments with either a coarse or fine nodalisation, and found that the flow instability in 
the experiments was not always predicted by the coarse nodalisation, and the source of the 
instability impacted whether finer nodalisation was required. 


Limitations of Constitutive Models 


Constitutive models are a key aspect of system codes, and are used to predict parameters such as 
interfacial heat transfer, interfacial drag, pressure losses, wall heat transfer and wall friction, based 
on the local flow conditions. The constitutive models are empirical or semi-empirical in nature and 
represent relevant interface transfer phenomena using spatial and time-averaged variables. While 
these models underpin the systems analysis approach, their limitations can have a significant 
impact on their ability to predict natural circulation (Roth and Aydogan, 2014a, Roth and Aydogan, 
2014b and D'Auria, 2017): 


« In general, each correlation has been developed to represent a single phenomenon using 
a limited set of experimental data (Volume 1, Section 4.6) usually under steady-state and 
fully-developed conditions, and therefore may have a limited range of applicability. Flows 
associated with natural circulation may deviate from these conditions, which may reduce the 
accuracy of the predictions and ability to predict the flow behaviour. 


The transition region between different flow regimes may not be well-defined, and each sys- 
tem code uses some form of interpolation technique. For typical forced circulation simula- 
tions, the impact will be small as frequent changes in the flow field are not seen. However, 
during natural circulation, flow oscillations can cause the flow field to change more frequently 
(and so the flow regime assumed will also change more frequently), and so the interpolation 
technique could have a larger impact. 


* System codes are not extensively validated for advanced reactor designs, such as non-water 
cooled reactors, and so the correlations may be less developed with a narrower range of 
applicability. 


For example, benchmarking work has indicated that RELAP5 may under-predict sub-cooled boiling 
at low power and low pressure (Shi et a/., 2018). This may affect the code’s ability to predict 
flashing instabilities in natural circulation, and make it difficult to demonstrate that no boiling occurs 
in the fuel channel under natural circulation. In this case, either subchannel codes could be used 
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to predict flow conditions in the core, or using conservative limits derived from subchannel and 
system code studies could be considered. 


In summary, code specific documentation should be reviewed to understand the correlations avail- 
able, any known limitations of the code, and whether any correlation is being used outside its range 
of applicability. The example above also highlights the benefit of literature review of similar cases, 
particularly for complex flows such as those that occur in natural circulation. 


Convergence 


System codes often feature automatic procedures to optimise the time step size to provide the best 
convergence and accuracy. The user usually provides minimum and maximum time step size and 
the code then determines the optimum value dynamically and adaptively. 


The solution methods employed tend to depend on the code, but most system codes employ a 
semi-implicit method and the choice of time step is limited by the Courant condition. In some cases 
(e.g. the nearly-implicit method in RELAP5-3D and 1D module in CATHARE), the need to satisfy 
the Courant condition is eliminated. Prior to advancing the hydrodynamic solution, the Courant 
limit check is carried out, and the time step size may be reduced. Due to the automatic adaptive 
selection of the time step size, solutions are in general insensitive to the user’s choice of maximum 
time step size. However, the automatic procedure may not always guarantee numerical stability. 


Time step sensitivity studies are recommended to ensure that the solution is not overly sensitive 
to the time step selection. Oscillations in the transient behaviour may be exaggerated by a finer 
discretisation. For example, in cases where flow reversal may be predicted within a time step, the 
user may need to enforce a smaller time step for accuracy. 


User Effects 


There are various user options when creating input files particularly for system codes (Petruzzi and 
D’Auria, 2008), including: 


* Physical model parameters (e.g. flow multipliers for choked flow, pressure loss coefficients 
and separator characteristics); 


¢ Nodalisation; 


* Specific system and components (e.g. scaled facility), specialised components such as pumps, 
valves and separators; 


* Initial conditions and boundary conditions. 


While this allows code flexibility to be applied to different reactor scenarios, it can result in large 
variations in results from different users when applying the same code for the same reactor and 
postulated accident scenario. These differences are considered as user effects, and are one of 
the most important sources of uncertainty. Sensitivity Analysis (SA) and Uncertainty Quantification 
(UQ) methods (as described in Volume 4) can be used to assess and quantify the uncertainty due 
to user effects. 


While the development of more advanced codes is supposed to reduce user effects, International 
Standard Problem (ISP) exercises initiated by the Committee on the Safety of Nuclear Installa- 
tions (CSNI) have shown that there are still significant user effects (D’Auria, 2017). Some ways to 
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minimise user effects have been proposed by CSNI groups, including user training and guidelines 
based on validation activities and quality assurance. 


Subchannel Analysis 


Subchannel codes are primarily used to predict the core flow behaviour and the occurrence of Crit- 
ical Heat Flux (CHF) in order to assess reactor operation against thermal hydraulic safety margins. 
The COBRA-TF family of codes are based on the two-fluid model with three fields (liquid, vapour 
and liquid droplets). The three-field model is shown to improve the code’s capability to predict en- 
trainment phenomena. However, other subchannel codes may employ different methods, such as 
the two-fluid model with two fields (i.e. liquid and vapour) and the three-equation approach (Cheng 
and Rao, 2015). 


Similar to system codes, some subchannel codes employ flow regime maps. For example, in 
COBRA-TF codes, the map is divided into a normal wall region and hot wall region: the normal wall 
region tends to correspond to pre-CHF region, while hot wall region corresponds to post-dryout re- 
gion. The flow regime in each mesh cell is determined using its fluid properties and flow conditions, 
and the appropriate models are selected to determine the closure terms including interfacial heat 
transfer, interfacial drag, wall heat transfer, wall friction and wall drag. 


Subchannel codes also allow the modelling of fuel assemblies such that full-core, pincell-resolved 
(i.e. one row of meshes/nodes per subchannel) simulations can be carried out. Whilst subchannel 
codes are, in principle, multi-dimensional there are significant differences between subchannel 
codes and conventional CFD analysis. These differences include: 


Computational requirements: In general, the number of cells used for subchannel analysis is 
of the order 10° to 10°. The number of cells used for CFD analysis is typically of the order 
10° so CFD is likely to be much more computationally expensive. 


Scale: Only component scale is available in subchannel codes while CFD can model different 
components and scales, if required. The minimum spatial resolution of subchannel codes is 
fixed by the size of a subchannel (usually order of 10~* m). 


Constitutive models: Subchannel codes are based on a number of empirical correlations, 
which reduces their accuracy (Section 3.1.2.2) and range of applicability. More complex flow 
behaviour can be modelled in CFD. 


Turbulence: Numerous turbulence models (e.g. RANS, LES, hybrid) are available in CFD. 
In subchannel codes, advanced turbulence models are not generally available. Since sub- 
channel codes assume that the flow is axially dominant, a simple turbulent diffusion model is 
normally used to compute turbulent transfer of axial momentum through channel gaps. 


As CFD is capable of analysing flow behaviour in fine detail (at increased computational expense), 
it becomes more appropriate for targeted analysis (e.g. analysing rod surface temperature for struc- 
tural analysis and pressure drop for spacer grids). 
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Nodalisation and Discretisation 


A subchannel code is often used to model domains, such as a fluid volume between rods, a lumped 
region of a core or a segment of downcomer. The subchannel can be viewed as a stack of mesh 
cells. The distribution of mesh cell length along the subchannel length can be uniform or non- 
uniform. If the subchannel includes a geometrical feature such as a spacer grid, a non-uniform 
mesh with smaller mesh cell sizes around the grid may be used (Avramova and Cuervo, 2013). 


Limitations of Constitutive Models 


Similar to system codes, some subchannel codes employ flow regime maps. Since subchannel 
codes are designed for vertical flow, it is common to exclude horizontal flow regime maps. 


Some codes allow users to choose between different models, such as the friction factor correlation 
(including Zigrang-Sylvester, Churchill and McAdams correlations). In this case, sensitivity anal- 
yses are recommended to determine an optimal choice of model. Some of the COBRA-TF code 
limitations (Salko et al/., 2016, Avramova et al., 2016) include: 


* COBRA-TF based codes tend to over-predict the rate of void generation and two-phase pres- 
sure. This may indicate that some improvements are needed in the interfacial, wall friction, 
sub-cooled heat transfer and turbulent mixing models. 


* Most codes experience difficulty in predicting void distribution near unheated regions. 


* Some codes have a bias with respect to pressure; there may be a tendency to over-predict 
critical power at lower pressure and under-predict at higher pressure. 


* Most V&V has been limited to single channel and small rod bundle configurations. 


¢ Subchannel codes alone are not equipped to capture all local effects of spacer grids (e.g. grid 
enhanced heat transfer, pressure loss, vane directed cross-flow and grid-enhanced turbulent 
mixing). Therefore, use of CFD (or coupling between CFD and subchannel codes) may be 
necessary if detailed modelling of the spacer grids is required. 


Convergence 


As for system codes, subchannel codes generally have automatic procedures to select the time 
step size for convergence and accuracy. The user provides the minimum and maximum time step 
size and the code will determine the time step dynamically and adaptively. 


In COBRA-TF based codes, there are two iteration loops per time step; the outer iteration for setting 
the continuity and energy equations over the entire mesh and the inner iteration for solving the 
pressure matrix which was created in the outer iteration. The maximum pressure change during 
the inner iteration is compared to the user-specified convergence criterion. If the convergence 
criterion is not satisfied, the code changes the time step size by taking into account the Courant 
limit, pressure change, void fraction change and error. It is noted that convergence issues can arise 
in low pressure conditions (< 4bar). 
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CFD Analysis 


This section considers modelling buoyancy affected flows using CFD. It builds on the more general 
discussion in Volume 1 (Section 4.5) and complements a number of more general sources of 
guidance for industrial CFD computations (e.g. industrial guidelines such as ERCOFTAC, 2000, 
NAFEMS, 2019, CSNI, 2015b and NAFEMS, 2003, and books such as Versteeg and Malalasekera, 
2007). 


As such, only aspects considered particularly relevant to passive cooling and natural convection 
are included. After introductory remarks, this section considers: 


* The CFD approaches available, aspects of mesh generation and case definition for these 
approaches (Sections 3.2.2, 3.2.3 and 3.2.4 respectively). 


* Specific aspects of using LES, RANS and hybrid approaches to model buoyancy affected 
flows (Sections 3.2.5, 3.2.6 and 3.2.7 respectively). 


« Judging convergence (Section 3.2.8). 


It is noted that while much of Section 2.2 is relevant to CFD analysis, Section 2.2.4 in particular dis- 
cusses key aspects of theory relevant to modelling buoyancy affected flows using CFD. In general, 
using CFD tools to solve natural convection flow fields is relatively complex, and best performed by 
experienced users. 


Introductory Remarks 


As noted in Volume 1, the motivation for the development of CFD methods extends widely beyond 
the nuclear industry, and CFD tools are well established in the design, analysis and assessment 
of engineering systems across a wide range of sectors. This greatly increases the pace of de- 
velopment of models relevant to challenging areas (such as buoyancy affected flow) and enables 
engineers to use the same tools and experience to address challenges in a range of industries. 


The computational cost associated with large-scale CFD means routine analysis of full passive 
cooling systems is likely to remain rare for at least the foreseeable future. System codes are de- 
signed for this and are often well validated for their intended use. However, buoyancy driven flows 
can exhibit complex behaviour even in simple geometries (Section 2.1), and closed passive cooling 
systems using natural circulation may exhibit unstable and transient behaviour (Section 2.4). 


In this context, CFD provides a means of predicting information that could not easily be obtained 
otherwise (particularly full field information for unsteady or complex flows) and can provide useful 
comparative data for system-level or experimental analysis (especially where local 3D flow features 
might be expected to significantly influence overall performance). CFD may have a key role in areas 
like detailed design and optimisation, exploring concepts, understanding detailed flow phenomena 
and supporting the planning and interpreting of experimental analysis. CFD is also used to provide 
information that can drive the development of lower fidelity models used in system codes and 
can help transfer experimental findings from reduced-scale to reactor conditions. As the cost of 
performing CFD simulations has fallen relative to the cost of experimental analysis, many industries 
are doing more CFD and less experimental work. 
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Whilst the complexity and applicability of CFD tools is vast, CFD predictions remain only as trust- 
worthy as the models they use and the user operating them. The models that CFD tools include 
for buoyancy effects may differ significantly in complexity and maturity, making the choice of model 
itself a complex task, requiring users competent both in the tool itself and in the underlying physics 
of the problem at hand. Volume 1 considers how CFD analysis might be planned, and notes that 
a staged approach to quality assurance is often beneficial, which is particularly true for technically 
complex work like the prediction of buoyancy driven flows. 


CFD Approaches 


A number of CFD approaches are available, primarily resulting from different ways of predicting 
turbulence, and these are introduced in Volume 1 (Section 4.5.3). The application of these ap- 
proaches to buoyancy affected flows is introduced briefly below and discussed in more detail in the 
following sections. There are a range of approaches and plenty of variations; for NTH, engineers 
must make informed judgements about the complexity and fidelity of their modelling work within 
the context of a graded approach informed by their safety and economic case. More general NTH 
guidance is provided in D’Auria (2017) and CSNI (2015b). 


DNS: Direct Numerical Simulation (DNS) can resolve all the length and timescales within a flow 
field. As discussed in Volume 1, it is considered highly accurate, but is very computationally expen- 
sive, so is predominantly used to study low Re or Gr flows in simple configurations. For buoyancy 
affected flows, DNS has been used to study basic scenarios such as fully-developed pipe flow 
with heat transfer and differentially heated cavities (Sebilleau et a/., 2018). The extremely detailed 
data that DNS can provide has then been used to develop fundamental understanding and refine 
the models used in other CFD approaches. This is particularly valuable for buoyancy driven flows, 
where experimental work is often challenging and there is often a lack of sufficiently detailed mea- 
sured data to use in CFD model development. However, in view of its limited industrial use, DNS is 
not considered further. 


LES: In LES, the Navier-Stokes equations are ‘filtered’, so that a proportion of the turbulence is 
resolved using the mesh and the remainder is modelled using a SGS model. As discussed in Vol- 
ume 1, LES is much less computationally expensive than DNS and potentially more accurate than 
RANS, but the computational costs (which scale with Re and Gr) are still high. Within the context 
of a graded approach for NTH, LES is most likely to be used in areas of high safety significance or 
commercial impact. For passive cooling systems, LES is most likely to be used to provide compar- 
ative data to gain confidence in RANS work, or to enable more detailed analysis of particular flow 
phenomena. 


RANS: In RANS, the instantaneous flow variables in the Navier-Stokes equations are decom- 
posed into mean and fluctuating parts, so that the mean flow is resolved and all turbulence is mod- 
elled (lower frequency unsteadiness in the mean flow can be captured using Unsteady Reynolds- 
Averaged Navier-Stokes (URANS), up to a point). As discussed in Volume 1, RANS approaches 
are used for the vast majority of industrial CFD (including for modelling buoyancy affected flows) 
and are often used to support more detailed analysis methods like LES or hybrid approaches. 
Within the context of a graded approach for NTH, RANS is therefore likely to be used in most sit- 
uations where CFD is needed, with the complexity of the modelling work reflecting the complexity 
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of the case and the safety significance and commercial impact associated with the work. Despite 
their well-known shortcomings for modelling buoyancy affected flows (which are discussed in the 
following sections and Section 2.2.4), RANS models can provide useful results with practical com- 
putational expense. For analysis of long plant transients containing buoyant flows, URANS is likely 
to be the most detailed CFD method that can be used with practical computational expense. 


Hybrid Methods: Hybrid methods combine RANS with aspects of LES within the same model, 
to take advantages of both. As discussed in Volume 1, there are a number of significantly different 
approaches available (and variants within these approaches). Hybrid methods may well be useful 
for modelling buoyancy affected flows, which often feature unsteady free-shear or highly separated 
flows away from walls. Resolving turbulent length scales in some areas of the domain is likely to 
increase the computational expense compared to RANS, but reduce it compared to LES. Like LES, 
the use of hybrid approaches is increasing in industry, and their use within a graded approach for 
NTH sits between the points made above for RANS and LES. 


Mesh Generation 


The majority of industrial CFD calculations use a computational mesh on which the flow is solved. 
The mesh determines the computational domain and is fitted to the geometry of interest. Its pur- 
pose is to enable gradients in the flow variables to be resolved. The mesh can affect the results of 
a calculation significantly (particularly for complex buoyancy driven flows) and may be difficult to 
update once work has started. 


Computational Domain and Geometry: These key aspects of CFD analysis are considered in 
Volume 1 (Section 4.5.2). Additional aspects that may be significant for buoyancy affected flows 
include: 


« Using symmetry planes or periodicity to reduce the size of the computational domain should 
be approached with care. Buoyancy affected flows are often asymmetric and aperiodic (even 
in symmetric or periodic geometry) as a result of the intrinsically unsteady and complex flow 
fields (such as plumes) that may occur. Using periodicity may also not allow enough space 
in the domain for flow effects to form. For example, in Rayleigh-Bénard convection between 
two large horizontal plates, the developing convection cells may have a diameter of the order 
of the domain height. However, if the domain width is insufficient, these convection cells will 
be unduly restricted and the predicted flow field may be incorrect. 


On simplifying geometry, studies have shown (IAEA, 2014, 2013a; Groetzbach, 2002) that 
natural convection flow fields can be particularly sensitive to thermal disturbances caused by 
structures and small geometrical features where flow is locally accelerated (such as the radii 
of hole edges). 
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Meshing Background: General background on mesh types and quality is provided in Volume 
1 (Section 4.5.2). More detailed guidance on creating meshes is provided in Volume 2 (Sec- 
tion 3.4.2). As noted in Volume 1, when a flow solution is available, the predicted flow should 
be visualised on top of the mesh and reviewed by an experienced CFD user who can consider 
whether the mesh is capable of appropriately resolving gradients that may (or perhaps should) 
exist in the flow. 


Natural or mixed convection flows may have complex features such as plumes or areas of recir- 
culation (Section 2.1) which will require sufficient mesh to resolve these regions of the domain. 
These will be in addition to other flow features which may be impacted by, or entirely unrelated 
to, buoyancy affects (such as boundary layers or areas of separated flow). Because of this, more 
mesh cells will generally be required in the interior of the domain than needed for forced convection 
flows, and mesh quality needs to be considered more carefully across the whole domain. 


If the mesh is too coarse, flow features may not be able to develop in the solution, leading to incor- 
rect flow predictions. For example, a plume has a shear layer around its circumference that may 
vary in space with time, which must be modelled to capture the flow feature. Ideally, enough mesh 
points are needed to model each resulting shear layer that may occur in the flow (typically at least 
10 normal to the shear layer, although software specific guidance should be checked). This may be 
challenging for industrial geometries, but can be investigated using validation data, sensitivity stud- 
ies and uncertainty quantification (Volume 1, Section 4.3). Specific aspects of meshing relevant to 
passive cooling systems are considered further in the following sections. 


Meshing Pools and Plena: CFD is used to predict flow temperatures in pools or storage build- 
ings, for example to check cooling flow is appropriately distributed. While pool geometries may be 
simple, the items stored within them may be geometrically complex. To avoid modelling these com- 
plex features, it may be appropriate to use a sub-model to predict the flow rates and temperatures 
associated with the heat input. Appropriate mesh is then needed both to model the heated flow 
entering the domain, and the resultant plumes. The resolution needed may be difficult to assess 
without prior knowledge, so initial calculations using periodic boundary conditions to roughly cap- 
ture part of the repeating geometry may be useful. A prismatic or hexahedral mesh is likely to be 
helpful within the volume to avoid dissipating plume flow features. 


For plena, the inflow and outflow pipework are often coupled with the plena to some extent, so it 
is likely to be worth extending the calculation domain down this pipework for some distance. Initial 
scoping calculations may be helpful to ensure that the boundary conditions are not over constrain- 
ing the predicted flow. The flow within plena themselves is often complex and three-dimensional, 
so providing sufficient mesh in the volume and capturing the flow transitioning from the pipework 
to the volume are likely to be important. 


Meshing Loops and Channels: For natural circulation applications, the flow in loops and chan- 
nels may be as complex as pools and plena, since the flow may not be aligned with the geometry 
and contain flow features like plumes. Meshes therefore are likely to be finer in the volume and 
have smaller aspect ratios than would be normal for forced convection, to enable buoyancy-related 
flow features to be predicted. 
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Case Definition 


Once a computational mesh has been developed, the CFD case can be configured by defining 
physical information and solver settings. This section focuses on the areas that may be particu- 
larly important for buoyancy affected flows and are general across CFD approaches (e.g. setting 
boundary conditions, fluid material properties and spatial and temporal discretisation). Turbulence 
and unsteadiness are also significant, and these are discussed separately in following sections. 
For buoyancy driven flows, the energy equation will need to be solved as well as the Navier-Stokes 
equations, because flow velocities and temperatures are intrinsically coupled. 


Boundary Conditions: Boundary conditions define the problem being solved, so if they contain 
inaccurate information or are not applied correctly, it is likely the results will also be inaccurate. 
The application of boundary conditions in CFD tools varies, so software specific documentation 
should be consulted. This is particularly true for buoyancy driven flows, where flows near boundary 
conditions can be more complex than for forced convection flows. For buoyancy affected flows, the 
following aspects are particularly significant: 


* Inflow and outflow boundary conditions should ideally be placed sufficiently far from areas of 
interest to avoid unwanted interactions (as noted in Section 3.2.3). If an inflow is fully devel- 
oped, suitable profiles for flow variables could be applied from literature sources, or created 
(e.g. by using a sub-model or extending the domain). Where a system is connected to ex- 
ternal or ambient conditions at both inlet and outlet (like the simple example in Figure 2.6), 
and those boundaries are represented by pressure conditions (rather than a specified ve- 
locity) then the solution may be harder to converge and care is needed to specify pressure 
correctly. The external hydrostatic variation between the two heights drives the flow, and this 
pressure difference should be applied. However, if the inlet and outlet are vertical planes, 
then a height varying pressure across the boundary would need to be specified. To avoid 
this, and to improve numerical convergence, CFD codes often allow an operating density 
to be specified, p,, and redefine the pressure field for the solution to be P’ = P — popgz, 
where z is the co-ordinate in the direction of gravity. In this case, a zero gauge pressure can 
be applied to all external boundaries and the hydrostatic contribution will be included. Care 
should be taken not to ‘double count’ by applying both an operating density and a prescribed 
pressure difference (this is a common error). 


¢ Inflow turbulence will need to be specified in an appropriate form for the modelling approach 
being used. In transitional-Re flow, such as may occur in natural convection, inflow turbulence 
levels can have a significant impact on results. For RANS, the values being specified can be 
converted into a turbulence intensity and length scale to check they are consistent with the 
case. For LES, the velocity specified at the inlet is necessarily unsteady and must include 
turbulent fluctuations. There are a number of methods to generate this ‘synthetic turbulence’ 
(most CFD codes provide at least one), but the resulting boundary conditions should be 
checked to ensure they have realistic statistics (see Jarrin et a/., 2006 and di Mare et al., 
2006 for further detail). 


« Wall boundary conditions may have a significant impact, because they are often used to 
introduce heat to buoyancy driven flows. For natural circulation loops, the fluid’s ability to 
exchange heat with solid walls, the thermal mass of these solid walls and changes in ambi- 
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ent temperature can affect flow stability (Basu et a/., 2007). In addition, even well-insulated 
systems may have important heat losses and axial conduction along pipework may be sig- 
nificant. As such, simply assuming adiabatic walls may result in incorrect behaviour. Either a 
uniform heat flux or uniform temperature boundary conditions are often applied, but consid- 
eration should be given to how appropriate this assumption is. For buoyancy affected flows 
consideration should be given to using a CHT approach that solves the solid temperatures 
at the same time as the flow, using a thin wall, shell conduction or fully meshed solid (Craft 
et al., 2010). This is discussed further in Volume 2. Near-wall treatments for RANS turbulence 
models are considered in Section 3.2.6. 


Fluid Material Properties: Fluid properties are particularly important for natural convection, as 
the flow is driven by changes in these properties. As such, sourcing and using properties is dis- 
cussed in some detail in Section 2.3. 


In order to simulate buoyancy effects, it is necessary to include gravity and variable density within 
a CFD model. For closed domains (i.e. there are no inflows or outflows), this needs careful consid- 
eration in order to ensure that the CFD solver conserves mass in the computational domain. If the 
average density in the domain is not calculated correctly each iteration, it can lead to errors in the 
mass conservation. Two common solutions to prevent this problem are: 


1. Providing the domain with at least one inlet/outlet. In a closed natural circulation loop for 
example, a constant pressure boundary might be added using a small bore T-junction along 
the top horizontal leg, effectively mimicking a physical thermal expansion vessel (which would 
be present in a real loop) without having a significant impact on the flow. This boundary then 
allows small amounts of fluid to enter and leave the domain as needed. 


2. Using the Boussinesq approximation for buoyancy. As discussed in Section 2.3.2, this ap- 
proximation is valid only for small density variations and thus, for most liquids, small temper- 
ature differences. For real passive cooling systems this may not hold, so the validity of this 
approach should be checked before use. 


Regardless of the solver, for steady-state situations, the Boussinesq approximation may well be 
useful in getting CFD calculations to start in a stable manner, before the approach is switched to 
another method once the flow field has developed from the initial guess. 


Spatial and Temporal Discretisation: Discretisation errors arise as a result of replacing the 
continuous derivatives in the governing equations with discrete numerical approximations that can 
be solved. Spatial discretisation provides a numerical approximation of the flow variable deriva- 
tives, while temporal discretisation provides a numerical approximation of the time derivatives in 
an unsteady (also called transient or time-resolving) calculation. The magnitude of discretisation 
errors depend on: 


* The nature of the flow. 

¢ The nature and order of the discretisation scheme. 

¢ For spatial discretisation, the size, quality and type of mesh employed (Section 3.2.3). 
* For temporal discretisation, the time step used (Section 3.2.5 and 3.2.6.2). 


Spatial discretisation errors generally manifest as an additional source of diffusion in the flow 
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(called ‘numerical diffusion’). This can potentially dominate over real molecular or turbulent dif- 
fusion, especially where meshes are coarse. Temporal discretisation errors can cause inaccurate 
prediction of unsteady flows, or prevent unsteadiness from occurring entirely. Mixing and unsteadi- 
ness are often significant areas of interest for natural convection. 


In general, low order schemes are likely to be more stable but cause more diffusion than higher- 
order schemes, and so are often used initially to stabilise and converge solutions (Section 3.2.8). 
For RANS, second-order accuracy is most often used (particularly for pressure and momentum 
fields) but first-order may be appropriate in some cases (e.g. turbulence fields). For LES, second- 
order or higher is often recommended (Benhamadouche, 2017). However, CFD solvers may offer a 
number of numerical schemes and specific documentation should be consulted for recommended 
settings. 


Sensitivity studies considering mesh resolution (and numerical schemes) can be used to assess 
the impact of spatial discretisation errors if required. Likewise, sensitivity studies considering the 
time step (and time discretisation scheme) can be used to assess the impact of temporal errors 
in transient calculations. Ideally, a combined spatial-temporal sensitivity study (considering spatial 
and temporal aspects together) should be used, but due to the large computational expense of this 
it is more usual to perform a time step sensitivity study using a mesh that has been determined to 
provide sufficient spatial resolution. Further guidance on understanding discretisation errors can 
be found in ERCOFTAC (2000) and CSNI (2015b). 


LES Aspects 


LES is introduced in Sections 3.2.2 and 2.2.4.1, meshing aspects are considered in Section 3.2.3, 
and this section presents some advice for conducting LES work. Further details on LES methods 
are available in Georgiadis et a/. (2010), Bouffanais (2010) and Meyers et al. (2008). Despite 
computational costs (discussed below) the application of LES to NTH is steadily increasing. Some 
industrial examples of LES include: 


¢ T-junctions (Smith et a/., 2011; Héhne, 2014). 

« PWR plena (Simoneau and Champigny, 2008). 

« Pressurised Thermal Shock (PTS) (Loginov ef a/., 2011). 

¢ Rod-bundles (Mikuz and Tiselj, 2016). 

¢ Natural convection in the primary vessel (Merzari et a/., 2017). 


Some further examples of modelling buoyancy affected flows in the academic literature include: 


* Differentially heated cavities (Ammour et a/., 2013; Kumar and Dewan, 2016). 

¢ Buoyancy-aided channel flows (Duan and He, 2017; Lau et a/., 2012). 

* Transition in natural convection (Padilla and Silveira-Neto, 2008). 

* Mixed convection around large spent fuel storage cylinders (Champigny ef a/., 2007). 
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Computational Cost: As noted in Section 3.2.2, the computational cost of LES is less than DNS 
but significantly higher than RANS. This is particularly so for high Re or Gr boundary layer flows, 
and because near-wall turbulence levels are coupled to surface heat transfer, the computational 
cost also increases as Ra increases. For natural circulation flows, the computational expense of 
LES is still likely to be prohibitively large for industrial calculations, limiting its use to studies of 
specific areas of interest. 


High Performance Computing (HPC) facilities are typically required to produce LES results within 
reasonable time frames. Even so, the computational demands of performing properly resolved LES 
will likely make it unsuitable for routine large scale industrial calculations for some time. Projections 
in Larsson and Wang (2014) suggest that LES is unlikely to replace RANS for design and optimi- 
sation within the next 30 years at least. Hybrid approaches, considered in Section 3.2.7, may offer 
a useful compromise between more accurate LES and less computationally expensive RANS. In 
particular, Wall Modeled Large Eddy Simulation (WMLES) is so commonly used to avoid resolving 
turbulence near walls it is often presented as an option for LES rather than a hybrid approach. 
However, for long plant transients, even URANS may be impractical in some cases (Section 3.2.6). 


Time Step: LES calculations are transient, so a time step is generally used by the solver. It is 
important that the temporal resolution should match or exceed the spatial resolution (i.e. the time 
step should be small enough to adequately resolve the flow as it passes through each cell). CFD 
software may incorporate a number of different tools to apply and vary the time step so specific 
documentation should be reviewed. 


The time step needed (At) may be estimated from the mesh size using the Courant-Friedrichs 
Lewy (CFL) number (CFL = UAt/Ax, where U is the flow speed and Ax is the cell length in the 
flow direction). CFL numbers of 1 or less are normally used for LES. Like the mesh size (Volume 
2, Section 3.4.2.4), the time step is often estimated using an initial RANS calculation. To ensure 
that CFL < 1 in each cell of the LES calculation, this is normally achieved by conservatively con- 
sidering a CFL = 0.5 (i.e. using At ~ Ax/2Upawns) to account for any differences between LES 
instantaneous and RANS averaged velocities, and any effects specific to the RANS solution. 


Care is required when starting an unsteady LES simulation from a steady RANS solution (or chang- 
ing a time step during a calculation) to ensure that a proper transient solution is achieved before 
starting to average the statistical data. It may take a long time for the solution to reach a statisti- 
cally averageable state (Section 3.2.8), and could require a large number of residence/through-flow 
times (usually between 3 to 10). The residence time can be estimated using the bulk flow average 
velocities or locally generated streamlines (from a point, line or plane). 


Sub-Grid Scale Models: A number of SGS models are available for LES, and their purpose 
and background are introduced in Section 2.2.4.1. As previously mentioned, although the different 
SGS models have their strengths and weaknesses, the prime benefit of LES arises from its ability 
to resolve, rather than model, a significant proportion of the turbulent spectrum. The main features 
of the flow are usually well represented and the SGS models are generally less complex and less 
significant to the accuracy of the simulation than RANS turbulence models. Some common SGS 
models are introduced in Section 2.2.4.1. To summarise: 


* The Smagorinsky model is regarded as simple and robust, but in view of its limitations the 
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Dynamic Smagorinsky model is generally preferred. However, while this does not require the 
Smagorinsky coefficient to be specified in advance, it inherits the weaknesses of the more 
basic model. 


¢« The WALE model is also straightforward, but provides better modelling near walls and in 
regions of laminar flow. As such, the WALE model is more likely to be used for the wall- 
bounded or low Re flows associated with natural circulation. 


Other SGS models may be available in specific CFD software, so software-specific documentation 
should be consulted for guidance. It is noted that WMLES and similar hybrid approaches (intro- 
duced in Volume 1 and discussed in Section 3.2.7) may be presented in a similar manner to an 
SGS model in software documentation, as these approaches are commonly used in industrial LES 
calculations to overcome onerous near-wall meshing requirements (discussed in Section 3.2.3), 
particularly where Re is not low. 


Transition: It is possible to have fully-turbulent, transitional and laminar regions present within 
the same buoyancy-driven flow field, and this may or may not have a significant impact on the 
accuracy of CFD predictions. There are also a number of different transition mechanisms (Sec- 
tion 2.2.3). As a result, predicting laminar-turbulent transition using CFD is important, but complex, 
and an active research area. LES can be used to predict transition, but is likely to be too expen- 
sive to be used for this purpose for most industrial applications for some time. While the cost of 
LES calculations is lower than DNS, the predicted flow behaviour is likely to be sensitive to the 
performance of the SGS model and the detailed specification of upstream turbulence (which can 
be difficult to obtain in industrial contexts). 


RANS Aspects 


RANS is introduced in Sections 3.2.2 and 2.2.4.2, and meshing aspects are considered in Sec- 
tion 3.2.3. This section presents some advice for conducting RANS modelling work. 


Steady and Unsteady RANS 


One of the key aspects of RANS is the ability to run models as either steady-state (RANS) or un- 
steady (URANS). While turbulent flows are, strictly, always unsteady (due to turbulent fluctuations 
in flow quantities), large-scale unsteady motion can also arise from non-turbulence phenomena! 
or be imposed externally by boundary conditions which vary in time. This large-scale ‘coherent’ un- 
steadiness is often also present in laminar flows and can be considered distinct from the unsteady 
small-scale ‘incoherent’ unsteadiness typically associated with turbulence. 


URANS methods exploit the above distinction by solving a form of the RANS equations which retain 
the time-derivative term. Large-scale low frequency unsteadiness is then resolved, and turbulence 
is modelled using the same models as steady RANS (for example, moving plume structures may 
be predicted that provide a mixing mechanism in addition to that provided by turbulence). This 
approach is well established and has been applied successfully to many flows involving buoyancy 
(Ammour et a/., 2013; Kenjeres and Hanjalié, 1999; Wilson et a/., 2015), and has been extensively 
reviewed (Benhamadouche, 2017; Kenjeres and Hanjali¢é, 2009; Speziale, 1998). 


Such as plumes, periodic vortex shedding behind a bluff body, flows subject to system rotation, or the convective roll cells 
that develop in Rayleigh-Bénard convection. 
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Steady-state calculations account for the majority of engineering analysis, but must be used with 
care for buoyancy affected flows because, like the plume rising from a cup of tea, natural and mixed 
convection flow fields are often unsteady even with steady boundary conditions (Section 2.4). Ob- 
taining a converged steady solution where flows are unsteady in reality may be possible (Sec- 
tion 3.2.8), but the predicted flow fields may not be meaningful, or depend on the number of itera- 
tions and applied under-relaxation parameters. Therefore, a steady RANS solution may not give 
the same results as a long-term time-averaged URANS solution. 


Moving from RANS to URANS can result in significantly different flow field predictions (Figure 3.1). 
If RANS and URANS predicted flow fields are the same, this may suggest that the real flow is 
broadly steady. 


a) Steady RANS b) Time-averaged URANS 


Figure 3.1: Example of impinging jet, RANS vs URANS (Contours of velocity magnitude 
(m/s) in TALL-3D test section from Study A). 


URANS provides a statistical representation of the flow as it evolves in time. It will not provide a 
time-exact reproduction of the flow that might be otherwise obtained with experimental, or even 
higher fidelity (LES, DNS) numerical, methods. In particular, caution should be used in interpreting 
coherent structures predicted in URANS (Benhamadouche, 2017). This should be considered if 
comparing URANS with LES or DNS. Some additional aspects of using URANS are: 


« Unlike LES, contributions from the turbulence model are not tied to the mesh-size and thus 
decreasing the mesh-sizing will not enable additional, turbulent, motions to be captured. This 
means that it is possible to obtain a truly mesh-independent solution using URANS methods. 


¢ Often, unsteadiness arises from instabilities in the flow which subsequently develop into co- 
herent motion. As turbulence is, however, principally a dissipative phenomenon, a turbulence 
model which is overly dissipative may prevent these instabilities from developing in the first 
place, suppressing unsteady motion. Generally, turbulence models which are less dissipative 
(typically those able to account for some level of anisotropy, so RSMs or non-linear EVMs) 
will be better at reproducing unsteady behaviour. This point reveals a subtlety with URANS; 
whilst the turbulence model is only responsible for modelling incoherent turbulent motion, a 
change in turbulence model may enable previously damped out unsteady coherent motion 
to be subsequently captured. 


* URANS assumes there is scale separation between the turbulent fluctuations and unsteady 
coherent motions’. As the frequency of the coherent motion increases, the distinction be- 


2 This implies an assumption that the temporal average of the turbulent quantities is not affected by the global unsteadiness. 
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tween what constitutes coherent motion and what constitutes turbulent fluctuations can be- 
come more difficult to discern. In such cases it may still be possible to obtain a converged so- 
lution, so care is required when interpreting results. In buoyancy driven flows, where plumes 
and large-scale roll cells may interact with both highly turbulent and laminar regions within 
the same flow, this distinction can be difficult to ascertain without analysis. 


As noted in Volume 1, for many decades URANS is likely to be the only practical CFD method 
for predicting long duration problems such as reactor shut-down transients (which may extend 
over several hours and include buoyancy affected flows). Given the points on unsteadiness noted 
above, URANS may well be the most appropriate method for solving natural convection flow fields. 
However, it is more complex than using RANS, and this should be considered when the approach 
to Verification, Validation and Uncertainty Quantification (VVUQ) is developed. Further details sur- 
rounding the meaning and interpretation of URANS methods can be found in Spalart (2000). 


Time Step 


If transient data is required from a URANS solution, a time step is normally specified®. The choice 
of time step will depend on the timescales present in the flow, and can have a significant impact on 
results. Identifying and estimating these timescales normally requires engineering judgement, but 
the below considerations may be helpful: 


* Explicit solvers generally need CFLmax ~ 1, while implicit solvers are usually less sensitive 
to numerical instability and so larger CFL values may be acceptable. Recommended values 
for CFL numbers are often provided in specific documentation for CFD software (and are 
particularly important for explicit methods). 


¢ For features or regions within the computational domain, a time step may be derived using 
At = t/N where N is an appropriate number of time steps to resolve motions on an esti- 
mated characteristic timescale (t). The number of time steps required is likely to depend on 
the numerical scheme (and number of iterations per physical time step for implicit schemes), 
so software-specific documentation should be consulted. Methods for estimating a charac- 
teristic timescale include: 

— For transient developing flows (e.g. a jet starting across a volume), it may be appropriate 
to use a characteristic flow velocity (e.g. the jet velocity) and the length over which the 

flow is expected to develop. 


— For buoyant flows (High Ra number), it may be appropriate to use a buoyancy-driven 
reference velocity t = L/U = L/.,/gBATL, where g is acceleration due to gravity, B is 
thermal expansion coefficient, AT is a characteristic temperature difference and L is a 
characteristic length scale consistent with this temperature difference. 


— For unsteady flows, it may be appropriate to estimate the expected flow oscillation fre- 
quency using the Strouhal number (Sr), and hence the timescales for vortex shedding. 


— For analysis supporting experimental work, it may be appropriate to estimate a time 
step based on experimental results. 


3 Some solvers may provide methods of obtaining time-averaged flow fields without a time step, but such results may not be 
helpful in understanding an unsteady buoyant flow field (for example switching behaviour between two states). 
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« Where aspects such as the boundary conditions, parameters or geometry are changing 
through time, a time step might be also be defined by At = t/N where N is a number 
of time steps needed to capture the given change occurring over a timescale t. 


Ideally, the chosen time step would be the smallest value predicted for the various aspects iden- 
tified. However, this may not always be possible due to the number of time steps that would be 
needed (e.g. if the smallest timescales are of the order of tenths of seconds but overall plant tran- 
sient timescales last for weeks) leading to a judgement about what timescale can be used. In some 
cases, it may be appropriate to use a system code to predict the plant transient behaviour, with a 
CFD calculation performed assuming conditions are quasi-steady-state at key points through the 
transient, or to train a correlation or sub-model using a number of CFD studies to build into the 
system code directly (Volume 2, Section 3.1). 


A time step sensitivity study should be performed to ensure that the chosen time step is appropriate 
for the case. Unlike LES or DNS, a URANS solution can, in principle, be independent of time 
step (as well as mesh-size). Starting an unsteady solution or changing the time step during a 
solution must be done with care, and like LES the time taken for the solution to reach a statistically 
averageable state should be considered (Section 3.2.8). Other considerations may be needed for 
solids in CHT cases (Volume 2, Section 3.4.5). 


Turbulence Models 


As introduced in Volume 1, turbulence modelling is a key aspect of the RANS approach. It is 
generally accepted that a ‘universal’ turbulence model ideal for all flows does not exist. However, 
turbulence modelling for single-phase CFD has matured to the point where reasonably reliable and 
good results can be obtained for a wide range of engineering applications with current computers. 


The choice of turbulence model should be made on a case by case basis, and is likely to involve 
aspects such as common practice, user experience, sensitivity studies, validation data, computa- 
tional resources, implementations within the CFD tools being used, as well as the technical aspects 
of the models themselves. A number of common models and their technical background are intro- 
duced in Section 2.2.4.3. For buoyancy affected flows, either a two-equation linear EVM (probably 
a variety of k -€ model or the k- w SST model) or a variety of RSM is likely to be chosen, although 
they are likely to be challenged by buoyancy affected flows. More details on using these models 
are provided in the following sections. 


The implementation of turbulence models in CFD tools can vary, so itis always necessary to consult 
software specific documentation (e.g. user/modelling/theory guides) for detailed information on the 
configuration and use of the modelling approaches offered. A number model variants are normally 
offered, and complexity in turbulence modelling is such that a variant that works well in some 
scenarios may fail in others in ways that can be difficult to predict. This is partly why VVUQ, which 
is introduced in Volume 1 (Section 4.3) and discussed in Volume 4, is a key part of industrial CFD. 


Turbulence models are empirically calibrated to a wide range of flows by the developer, which may 
or may not be directly applicable to a given case. However, recalibrating turbulence model coef- 
ficients needs great care because it is possible to void the original calibration, cause unintended 
or unphysical behaviour, mask unrelated modelling deficiencies and lose the ability to compare 


54 of 91 


Methodology 


results with other work. Such recalibration would need detailed validation (e.g. using detailed ex- 
periments and/or DNS) and may therefore be very expensive. However, some turbulence models 
(such as the GEKO model, Section 2.2.4.3) provide additional parameters that enable tuning within 
controlled bounds. While this may offer benefits in some flow regimes, these parameters should be 
used with caution as these models are fairly new and have little published application to buoyancy 
affected flows. 


In summary, particularly for buoyancy affected flows, the turbulence modelling approach chosen 
will be dependent on the individual case, software-specific documents should be consulted, and 
sensitivity to different modelling approaches may be tested as part of VVUQ. 


k-e and k-w: These two-equation linear EVMs are by far the most commonly used models 
in industrial CFD, and are widely used in NTH applications. They are simple to use and com- 
putationally affordable. They may be particularly challenged by natural or mixed convection flows 
(Section 2.2.4.3), but given their well-known limitations it can be surprising how well they perform 
in complex industrial cases, including mixed and natural convection. Global quantities, such as 
pressure drops, mass flow rates, bulk temperatures and averaged Nusselt numbers etc., can be 
reliable (Benhamadouche, 2017), although some turbulence quantities are likely to be less accu- 
rate. Whether this is acceptable or not will depend largely on the application. Where more complex 
models are used (such as RSM or LES), these models are often helpful to gain initial solutions. 
More advanced EVMs have been proposed (Section 2.2.4.3) but are little used. Example uses of 
two-equation linear EVMs include: 


Vertical mixed convection flows (Cotton and Jackson, 1990; Keshmiri et a/., 2012, 2008). 
Natural convection in enclosures (Altag and Ugurlubilek, 2016; Hsieh and Lien, 2004; Mirosh- 
nichenko and Sheremet, 2018). 

Differentially heated cavities (Omranian et al., 2014; Rathore and Das, 2016). 


Natural circulation loops (Naphade et a/., 2013; Wang et a/., 2013; Kudariyawar et a/., 2016). 
Ribbed passages with heat transfer (Keshmiri et a/., 2016). 

Lower core of a PWR (Fournier et a/., 2007). 

Reactor cavity of modular HTGR (Zhao et al/., 2017). 

Passive residual heat removal heat exchanger (Ge et al., 2018; Zhang et al., 2015). 
Internally heated melt pools (Mao et al., 2018). 

Heat transfer due to impinging jets (Craft et a/., 1993; Sharif and Mothe, 2009). 

Buoyant plumes (Kumar and Dewan, 2014; Macpherson and Tunstall, 2020). 


k-e: It can be difficult to choose between the three main variants in this family of models. In 
principle, RNG or Realizable should be at least as accurate as the standard k - € model for many 
applications. However, mixed results have been reported in the literature for buoyancy affected 
flows and thus neither variant provides a clear and consistent advantage over the standard or low- 
Re, form of the k-¢ model (Chen, 1995; Faheem et a/., 2016; Grassi and Testi, 2007). In fact, the 
low-Re form of the standard model (Launder and Sharma, 1974) has demonstrated reasonably 
good performance in predicting heat transfer in several buoyancy influenced flows, including verti- 
cal mixed convection (Keshmiri et a/., 2008) and differentially heated cavities (Ince and Launder, 
1989; Omranian et al., 2014). Further comparisons of various low Reynolds k -< models can be 
found in Costa et al. (1999). 
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k-w: The k-w SST model is the most popular variant of this family of models. It has proved 
successful in a number of situations, including external aerodynamic and internal wall-bounded 
flows with adverse pressure gradients. Less attention has been paid to understanding its use and 
reliability in buoyancy affected flows compared with k-e models, which might arguably reduce 
confidence, although good results have been reported for wall-bounded plumes (Macpherson and 
Tunstall, 2020). 


RSM: Despite providing a more sound theoretical basis for modelling than EVMs (Section 2.2.4.3), 
RSMs have not been as widely used in industrial CFD as once anticipated. A variety of factors 
are thought to contribute towards this, including the increased number of coupled equations re- 
ducing numerical stability and increasing computational cost. Some early studies may have been 
compromised by low mesh resolution making it hard to differentiate between model performance 
and numerical errors (Benhamadouche, 2017). However, a number of more recent studies have 
demonstrated clear improvements: 


¢ PWR primary loops: Bellet and Benhamadouche (2011) showed that even advanced EVMs 
(such as the v?- f model) predicted the tangential velocity profile in a confined vortex tube 
poorly compared with the SSG RSM. 


« Inclined heated cavities: In an unstable configuration, Omranian et al. (2014) demonstrated 
many models, including EVMs, could return reasonable results since the flow is dominated 
by large-scale unsteady roll cells. However, in a stable configuration, a RSM together with 
advanced wall functions was needed to capture the subtle secondary flows present in the 
measured flow field. 


« Buoyancy driven counter-current flow within a pipe: EBRSM was used by Sebilleau ef al. 
(2016) to predict buoyancy driven counter-current flow, demonstrating improvements over 
both k-¢ and k-w linear EVMs. 


« PWR plenum: Martinez and Alvarez (2009) demonstrated that secondary motions are well- 
captured by a RSM once a fine enough grid is employed to adequately represent the intricate 
geometric details. 


Since RSMs may be less numerically stable than EVMs (especially when starting from relatively 
unrealistic initial conditions), it is often useful to start the calculation with an EVM before switching 
to a RSM when the predicted flow is more realistic. RSMs are also more sensitive to discretisation 
errors than EVMs (Benhamadouche, 2015), so a fine mesh and at least second-order discretisation 
scheme may be needed (see Section 3.2.4). First-order schemes may, however, improve numerical 
stability at the start of the calculation. 


Code documentation should be checked to ensure buoyancy effects are properly included, as this 
can be inconsistent’. In addition, not all RSMs are valid within the near-wall region, so documenta- 
tion should be checked to ensure that models are compatible with the near-wall modelling approach 
chosen. 


4 In one commercial code, buoyancy effects are included with an e-based RSM but cannot be included with an w-based RSM. 
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Near Wall Modelling 


For natural convection, the modelling of the near-wall region can have a significant impact on model 
predictions. As introduced in Volume 1 (Section 4.5.3), RANS near-wall modelling approaches 
typically include wall resolving (low-Re, fine mesh, small y*), standard wall functions (high-Re, 
coarser mesh, larger y*), enhanced wall functions (CFD code tailors approach to local y+) and 
advanced wall functions (more complex models catering for non-equilibrium conditions), although 
the terms used and detailed implementations may vary between CFD software. 


A wall resolving approach is likely to provide more accurate solutions than a standard wall function 
approach for buoyant flows (see below). In addition, natural convection usually results in much 
lower flow rates than in forced convection, and thus the computational cost of a wall resolving 
approach may be reduced. However, in many industrial contexts, obtaining a mesh to provide a 
small enough yt across the whole domain for a wall resolving approach is likely to be difficult 
and lead to unnecessary computational expense. Therefore, an enhanced wall function approach 
is often used in industry, with the mesh designed to enable near-wall flows to be resolved in key 
areas and wall functions used elsewhere. Developing a mesh to enable this is considered in Volume 
2 (Section 3.4.2), which considers near-wall meshing generally. 


If a wall-resolving approach is used (either on its own or through enhanced wall functions), the 
turbulence model must be valid for use in the near-wall region. w-based models are valid near 
walls, but low-Re versions that improve near-wall predictions may be available. Many original «- 
based models are not valid near walls, but many low-Re extensions exist (such as Launder and 
Sharma, 1974) and are widely used. Some low-Re turbulence models may use functions based 
on wall distance, which may be a problem in complex geometries (wall distances may be hard to 
define). Other options (e.g. based on turbulent Re) could be considered if available. 


Standard wall functions based on the logarithmic law-of-the-wall should be used with care, as they 
were developed for fully-developed flows subjected to simple shear, with the near-wall region in 
local equilibrium. This is highly unlikely to be true in buoyancy-driven boundary layers, where turbu- 
lence creation and destruction rates are far from balanced and transport effects may be substantial. 
These limitations are also shared by the corresponding formula for the thermal field. Standard wall 
functions may therefore cause significant errors in natural or mixed convection flows, especially 
where transition might be expected. 


A number of advanced wall functions have been developed to account for buoyancy effects (see, for 
example Craft et a/., 2002), which have shown good performance in natural and mixed convection 
flows (Craft et a/., 2006; Omranian et a/., 2014). Extensions to account for the effects of rough walls 
have also been published (Suga et al., 2006). However, these are not currently widely available in 
commercial CFD solvers. 


Turbulent Heat Transfer Models 


The modelling of turbulent heat fluxes can have a large impact on predictions of buoyant turbulent 
flows because it can directly affect surface heat fluxes and mean temperature profiles®. The choice 
of model will largely depend on the approach used to model the Reynolds stresses, since these 


5 This is analogous to the impact of Reynolds stresses on the mean velocity profile. 
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appear directly in most advanced turbulent heat flux models®. CFD software codes may have lim- 
ited models/options available, and so the code documentation should be checked. A number of 
models are introduced in Section 2.2.4.4 and their use is discussed below: 


* If a linear EVM is used for the Reynolds stresses, there is generally little benefit in using 
more advanced models than the SGDH, although a modified version of the GGDH has been 
used with success in natural convection (Ince and Launder, 1989). SGDH is best suited to 
simple forced or mixed convection flows (where only the cross-stream component of the 
turbulent heat fluxes is significant) and is therefore likely to struggle in more complex natural 
convection flows (Kenjeres and Hanjali¢, 1999). 


* The GGDH is likely to be the simplest approach used with a turbulence model capable of 
reliably predicting individual Reynolds stresses (i.e. a non-linear EVM or RSM). Even then, 
its implementation in CFD codes may be simplified to improve numerical robustness. 


* DFMs will normally only be used with RSMs and for specialist work due to the large increase 
in cost and numerical stiffness (four or five extra transport equations in 3D). 


Further guidance is provided in Hanjali¢é (2002) and Kays (1994). 


Transition 


Turbulent to Laminar: High-Re RANS models are calibrated for fully-turbulent flows and thus 
usually cannot take account of the viscous effects which lead to laminarisation. However, the mod- 
ifications included in some low-Re models mean that these are capable of (and often quite suc- 
cessful in) capturing laminarisation. For the differentially heated cavity discussed in Section 2.2.3, 
low-Re k- € models may well predict laminarisation as a result of these modifications, although the 
physical mechanisms are not modelled. 


Laminar to Turbulent: This poses a significant challenge for many RANS based methods. Whilst 
they cannot inherently predict natural transition (since the instabilities which cause it are averaged 
out), some low-Re models have shown limited success in modelling bypass or separated transi- 
tion (Menter et al., 2006). Transition specific models (such as k-w SST +) introduce additional 
variables (i.e. turbulence intermittency) to improve the accuracy of transition prediction. However, 
these models are generally developed for predicting transition in classic boundary layers not in 
other wall-bounded flows, free-shear flows or fully developed pipe flows. Further, they may only be 
designed to predict certain transition mechanisms (e.g. bypass transition, Section 2.2.3) and may 
not be calibrated to predict transition in buoyant flows. 


These are clearly significant limitations for predicting transition in natural convection flow fields. 
These models should be approached with caution, particularly if transition is not clearly expected 
based on local geometry or flow features, or good quality validation data from a well understood 
system is not available. The existence of transition within a flow should be anticipated through con- 
sideration of the governing non-dimensional parameters and examination of any relevant experi- 
mental data. Approaches which use standard wall functions will not be able to accurately predict 
laminar-turbulent transition, and a low-Re approach is recommended as a minimum. 


More advanced models are unlikely to provide a benefit if the turbulence model cannot provide appropriately reliable 
Reynolds stresses. 
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A RANS/URANS simulation that has become fully laminar because of a lack of turbulence produc- 
tion (by mean strain or otherwise), cannot subsequently become turbulent, even if those production 
mechanisms do later become significant’. However, if only part of the flow laminarises, turbulence 
may be convected from other parts of the domain if the flow field allows. As a result of this, if the 
velocity field must be initialised at zero when modelling a buoyancy driven flow, the fields repre- 
senting the turbulence (e.g. k, €, w) should be initialised at such levels that turbulence does not 
completely decay before buoyant motions (i.e. mechanisms that generate turbulence) develop and 
subsequently sustain the turbulence. 


Hybrid Approaches 


There are many hybrid approaches available with various names that may require different meshes 
with different strengths and weaknesses, settings to use and subtleties in results interpretation. 
Some common approaches are introduced in Volume 1, including WMLES, Zonal or Embedded 
LES, the many variants of Detached Eddy Simulation (DES) and Scale-Adaptive Simulation (SAS) 
(Figure 3.2). There is no ‘ideal’ approach for all situations, and the best approach to use is likely to 
be application and software specific. It is therefore necessary to obtain a detailed understanding of 
the approaches available within a given software before using them. 


LES Zone RANS Zone 


RANS Zone 


LES Zone 


(a) Wall Modelled LES (WMLES) (b) Zonal LES 


RANS Zone RANS Zone 


— OF 
LES Zone Ae 


LES Zone 


(c) Detached Eddy Simulation (d) Scale Adaptive Simulation 
Figure 3.2: The main hybrid methods. 


Broadly however, some areas of the computational domain are likely to be solved using a turbu- 
lence resolving approach (often LES) and remaining areas solved using RANS, with flow informa- 
tion passed between the two. As such, much of the discussion provided in previous sections can 
be considered for the relevant part of the domain and flow solution. Specific to hybrid approaches: 


* It is necessary to understand where in the domain the scale resolving method is needed, 
so that the mesh can be tailored appropriately. This can be particularly difficult to assess 
for buoyancy affected flows without a flow solution, so an initial RANS solution may help 
identify areas to use scale resolving methods and support mesh sizing (Section 3.2.3) and 

7 Because the mathematical form of the production term, which is exact, contains the Reynolds-stresses themselves. Thus 


even if the mean velocity gradients later become significant they will be multiplied by the current (zero) value of the Reynolds 
stresses. Terms representing buoyancy generation are similarly affected, since they contain the turbulent heat fluxes. 
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selecting an appropriate time step (Section 3.2.5). Tailoring the mesh and time step to the 
flow is particularly important for LES approaches, such as WMLES, Zonal LES or DES. SAS 
may offer more flexibility on mesh sizing than more traditional LES-based methods, depend- 
ing on the level of detailed flow information required. As the details are model dependent, 
software-specific advice should be consulted. When a flow solution is available, this can be 
interrogated to check that the solver is being used correctly and consistently with the mesh. 


A key feature of many hybrid approaches is that information is passed between scale- 
resolving methods (often LES-like) and non-scale resolving methods (often RANS-like). This 
is particularly significant for turbulence information, and is often a key difference between 
approaches (particularly for DES variants) so it is necessary to understand how this is man- 
aged by consulting software-specific advice. It may be possible to apply synthetic turbulence 
at the inflow to LES aspects of the flow. This should be considered if the inflow is turbulent 
and it is appropriate to the models being used (consult software-specific advice). If used, it 
is normally necessary for the mesh resolution to be of LES quality to ensure this turbulence 
is transported properly. When a flow solution is available the performance of the solver and 
mesh around these interfaces can be considered. 


WMLES is commonly used in industrial LES to reduce mesh requirements near walls (Sec- 
tion 3.2.3), and is often necessary unless Fe is low. It may be possible to use this for the 
scale-resolving parts of other hybrid approaches as well; some approaches may need to use 
WMLES to obtain correct behaviour near walls (consult software-specific advice). 


Significant development of hybrid approaches is continuing, partly driven by large-scale projects 
(such as the European Union DESider project, Haase et al., 2009). Gritskevich et a/. (2014) ap- 
plied a number of different approaches to the OECD/NEA Vattenfall T-junction benchmark test case 
(Smith et al., 2011), including SAS, Delayed Detached Eddy Simulation (DDES) and Embedded 
LES (ELES). These models were able to accurately predict velocity profiles, although SAS results 
were more sensitive to settings than for DDES and ELES. However, the application of these meth- 
ods to buoyancy driven flows is much more limited (Kenjeres and Hanjalié, 2006; Ding et a/., 2019; 
Kocutar et a/., 2015) so there is less guidance available. 


While hybrid approaches may offer improved predictions, they are still relatively new and thus do 
not currently have the benefit of decades of extensive use associated with RANS (or pure LES). 
It is therefore appropriate to be cautious when using these methods, especially for flows strongly 
affected by buoyancy. Software-specific documentation (such as Menter, 2015) should be reviewed 
to gain a good understanding of the limitations associated with the specific model chosen. Further 
general advice can be found in Fréhlich and von Terzi (2008), D’Auria (2017) and a review by 
Holgate et al. (2019) provides details on recent advances. 


Convergence 


CFD calculations are generally iterated until the predicted flow field is considered converged (i.e. 
further iterations would give little benefit). There is no universal way of assessing convergence, 
and it may be particularly difficult to assess for complex flows affected by buoyancy. As such, 
experienced engineering judgement is often important, based on appropriate metrics. The metrics 
chosen are usually case-specific, but will typically include monitoring the detailed flow variables 
and consideration of residuals as well reviewing the predicted flow field itself (including checking 
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that the mesh is appropriate to resolve gradients in flow variables, Section 3.2.3). For unsteady 
calculations, locations where significant unsteady flow occurs are often visualised over a number 
of time steps to check that the time step is sufficient. 


Monitors: Relevant flow variables should generally be monitored at key locations in the flow. 
These point monitors should be placed using engineering judgement (e.g. in locations considered 
particularly important, or where strong variation in flow variables is expected). Integral monitors 
are also used to monitor variables such as mass flows or average temperatures across a surface. 
Both the value and behaviour of the monitored quantities should be considered in the light of 
engineering judgement and intuitive insight on expected behaviours. For point monitors, this might 
include comparing the steady-state values with expected behaviour, studying the nature of any 
oscillations and calculations of moving averages. For integral monitors, this might include checking 
mass and energy conservation through boundary conditions and within the domain. 


For steady-state calculations, the solution should be advanced until values no longer change with 
further iterations. It is noted that unchanging variables at monitor locations do not provide any 
indication of solution accuracy or that the solution has necessarily converged, only that further 
iteration may not change the monitor value (situations where this may be significant include CHT 
cases, or where under-relaxation is used heavily). In a complex case, it is likely that a large number 
of monitors will be needed (different flow variables at different points in the flow). 


For unsteady calculations, the simulation must be solved for a sufficient number of time steps 
in order to reach a statistically averageable state before generating any required statistical data. 
However, areas with relatively slow moving flow (such as are likely to occur in passive cooling 
systems) may take 3 to 10 residence times through the domain to reach an acceptable solution. 
This should be confirmed using monitors that should be located where key quantities are expected 
to vary the most, and may also support judgements on whether the time step used is suitable. 
For scale-resolving simulations (such as hybrid methods or LES), the mean flow parameters may 
converge more quickly than fluctuating parameters (such as the Reynolds stresses), which should 
be considered when choosing monitor variables. 


Residuals: Residuals are often scaled or normalised by physical quantities associated with the 
flow to make them more physically meaningful. Most CFD software does this automatically, but 
specific documentation should be consulted to ascertain their meaning and presentation before 
using them to judge convergence. In general, a reduction of at least three or four decimal orders of 
magnitude is often sought, although this may depend on how good the initial guess of the flow field 
used to initialise the calculation is (a good initial guess may be associated with smaller residual 
reductions, while a poor one may be associated with larger residual reductions). The behaviour 
of the residuals throughout the solution process should be monitored to ensure that they do not 
start to increase, which may indicate that the calculation is diverging. While residuals are easy to 
obtain and review, it is very unlikely that simply checking them will enable an adequate judgement 
on convergence, particularly for buoyancy affected flows. Flow quantities (local, integral or most 
likely both) should always be monitored in tandem with residuals. 
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Convergence Problems: Fixing convergence problems is normally case specific. The case should 
be checked (if not done already) to ensure it has been configured appropriately and consistently. 
The relationship between the predicted flow field and local aspects of the case such as the ge- 
ometry, mesh and boundary conditions can be investigated, and many CFD solvers enable the 
location of maximum residual to be plotted to aid this. Often problems can result from the mesh 
(perhaps as a result of the relationship between local mesh quality and flow variable gradients) 
or from inconsistencies in boundary conditions (usually between inlet and outlet). For steady-state 
RANS calculations, difficulty in obtaining convergence may indicate that the real flow is unsteady®. 
Switching to a URANS approach may allow the numerical method to capture unsteadiness in the 
mean flow variables and improve convergence (Section 3.2.6.1). For unsteady calculations, using 
a smaller time step (or solving more iterations per time step in iterative schemes) may be useful. 


In addition, first-order schemes and under-relaxation parameters within the CFD software can be 
used to enable the flow to converge initially, by reducing them from their default values. However, 
once the flow has stabilised, the flow should ideally be converged using second-order schemes 
with the under-relaxation parameters set to their default values. This is particularly true for the 
energy equation as the temperatures may still be changing (just extremely slowly), even though 
the residuals/monitors appear converged. 


Experimental Methods 


The role of experimental methods in NTH analysis and the value of experimental and modelling 
teams working together closely is introduced in Volume 1 (Section 4). This section builds on this, 
discussing flow rate and velocity measurements that might be used in passive cooling experiments, 
while temperature measurements are discussed in Volume 2 (Section 3.5). Flow visualisation and 
two-phase flow measurement techniques are discussed in Frazer-Nash (2019). 


Scaling 


As noted in Volume 1 (Section 4.6.5), it is rarely possible to perform experimental analysis of test 
rigs built to full scale, so some scaling is needed to ensure that the phenomena and performance 
seen in the experimental rig is representative of the full scale plant. This scaling may be challenging 
for passive cooling systems, as a result of the potentially large number of parameters that influence 
the performance of the system. A detailed introduction to scaling is provided by CSNI (2017). 


A detailed example of scaling for natural circulation is provided by Reyes and Hochreiter (1998), 
using a Hierarchical 2-Tiered Scaling (H2TS) approach to design an ‘ideally scaled’ Integral Effect 
Test (IET) facility to study natural convection. The Advanced Plant EXperiment (APEX) test facility 
was designed on the basis of this scaling analysis to provide a geometric representation of the 
Westinghouse AP600 nuclear steam supply system. 


The H2TS scaling approach was developed by Zuber (1991) to provide a comprehensive scaling 
approach that reduced the reliance on judgement to define scaling requirements, and consists 
of four stages: system breakdown, scale identification, top-down scaling analysis and bottom-up 


Especially if residuals, or monitors, appear to oscillate, although the frequency content of such oscillations may be mean- 
ingless as the solver is not necessarily time-accurate. 
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scaling analysis. It was intended to be widely applicable to engineering tests for the purposes of 
code validation and aligns the test priorities with the complex applicability structure of analysis 
codes. 


H2TS has been widely adopted for the design of more recent IET facilities and has been considered 
successful, primarily due to the comprehensive and auditable nature of the approach. As scaling 
distortions are inevitable for complex systems, the quantification of uncertainties resulting from 
scaling effects is a key consideration within the United States Nuclear Regulatory Commission 
(US NRC) Evaluation Methodology Development and Application Process (EMDAP) and as part of 
a Best Estimate Plus Uncertainty (BEPU) approach. 


Another example of scaling for natural circulation is provided by Ishii (2016), using a three-level 
scaling approach to design an ‘ideally scaled’ IET facility to study natural convection. This ap- 
proach considers firstly a global scaling analysis based on dimensionless numbers; secondly a 
scaling analysis based on the main system components and the interfaces between them; and 
thirdly analysis based on local phenomena. The ‘engineering scaled’ facility is then designed to 
approach the ‘ideally scaled’ facility as closely as possible, taking into account practical engineer- 
ing constraints (such as standard sizes of pipework and safety-related pressure limits that might 
apply to test environments). System codes are also used to predict the overall performance in 
both facilities under steady-state conditions, to improve confidence in the overall scaling and to 
demonstrate matching with other facilities considered representative of a real PWR design. 


While developing detailed scaling to provide performance data may be complex, the overall qualita- 
tive behaviour of the flow in systems may not be greatly affected by the scale of the testing facility. 
For example, Schultz et a/. (1987) reviews work on a number of PWR natural circulation rigs and 
observes similar overall system performance and detailed steam generator flows, despite the tests 
being at a large range of scales. 


Integral Measurements 


Integral measurements include the mass flow rate, temperature differences across any heat ex- 
changers and also spot temperature measurements at key locations within the system. These 
measurements are often used to understand the system performance as a whole and provide data 
to compare with system analysis (Section 3.1). 


Two key measurements required for assessing the performance of a passive cooling system are 
the mass flow rate and heat transfer rates of the heat exchangers (i.e. heater or cooler): 


¢ Sensitivity of natural circulation flows to disturbances mean that instrumentation should be 
as non-invasive as possible. This needs particular consideration for measuring the flow rate 
(Section 3.3.2.1), since conventional flow meters or orifice plates are likely to have a signifi- 
cant impact on the flow. Low impact methods include Doppler shift ultrasonic or electromag- 
netic flow meters, alternatively the flow rate can be deduced by measuring the pressure drop 
across particular sections of the loop using pressure transducers (e.g. over a horizontal leg 
or between the inlet and outlet of a heat exchanger). 


* For heat transfer measurements, both the heater and cooler heat flux should be quantified if 
possible. Probes should be positioned to provide information about both local behaviour (e.g. 
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inlet and outlet of heater) and the global system behaviour (at enough positions to provide a 
reasonable measure of, for example, the loop fluid temperature). 


Flow Rate Measurement Techniques 


A number of common measurement methods are summarised below. It is worth noting that many 
flow measurement techniques, particularly non-intrusive methods, rely on an assumption that the 
flow is fully developed. This normally requires at least an L/D > 10 from an inlet or bend, and may 
need an L/D = 100 for the flow to be really fully developed (Klein, 1981). This may not be met in 
a real plant situation, increasing the uncertainty of the measurements. 


Orifice Plate: An orifice plate is a thin metal disc containing one or more holes, which is inserted 
into the pipe carrying the flow (LaNasa and Upp, 2014). Flow is restricted by the hole which in- 
troduces a differential pressure across the orifice. The upstream and downstream pressures are 
measured using a differential pressure gauge, and a calibrated pressure loss coefficient is applied 
to calculate the flow rate. Orifice plates are a simple, cheap and widely used method for measuring 


volumetric flow rate. However, measurement inaccuracy is typically + 2%, but can reach as high 
as + 5% (LaNasa and Upp, 2014). 


Coriolis Flow Meters: These are typically used to measure the mass flow rate of liquids, but 
can also be used on some gases and can be applied to indirectly measure volumetric flow rate. 
A Coriolis flow meter contains a pair of parallel vibrating tubes, or a single vibrating tube formed 
into two parallel sections (Morris and Langari, 2020). As fluid flows through the vibrating tube(s), 
a Coriolis force is generated and the tube is deflected further to the existing vibratory motion. This 
deflection is proportional to the mass flow rate of the fluid, which can be measured using a suitable 


sensor. Coriolis flow meters may have a high accuracy of + 0.2% (Morris and Langari, 2020) and 
can be used on liquids, gases, slurries and two-phase flows (LaNasa and Upp, 2014). However, 
they are expensive and can suffer from mechanical problems such as fatigue and corrosion. 


Doppler Shift Ultrasonic Flow Meter: A Doppler shift flow meter consists of an ultrasonic trans- 
mitter and receiver clamped to the outside of a pipe or fluid carrying vessel (Morris and Langari, 
2020). The transmitter emits ultrasonic waves which are deflected by scattering elements in the 
fluid and received by the receiver. This deflection causes a change in wave frequency, which is 
used to determine the flow rate. Doppler shift flow meters are typically inexpensive and can be 
used on gases or liquids; however the measurement accuracy depends on a number of different 
parameters, and so accurate measurements require careful calibration. These instruments are typ- 
ically used for flow indications instead of accurate quantification of volumetric flow rate (Morris and 
Langari, 2020). Variants with higher measurement accuracy have been developed, but these can 
be significantly more expensive. However, since these flow meters are clamped to the outside of a 
pipe, no contact is required between the flow meter and fluid. 
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Electromagnetic Flow Meter: Electromagnetic flow meters use Faraday’s Law to determine the 
flow of liquid in a pipe, as a conductive liquid flowing through a magnetic field generates a volt- 
age. Since the voltage generated is proportional to the velocity of the liquid, the flow rate can be 
measured from the voltage signal. The main advantages are that it is non-intrusive with no moving 


parts and a typical accuracy of up to + 0.5%. However, they can be relatively heavy and only work 
for conductive fluids, and so gases and deionised water cannot be measured. 


Detailed Flow Field Measurements 


Passive cooling systems often use natural convection phenomena to drive the flow and thus de- 
tailed measurements of such phenomena in isolation, through use of idealised loops or very simple 
geometries (differentially heated cavities, for example), can advance both physical understanding 
and the development of numerical models. Models validated in these simple scenarios can then 
be taken forward into more complex situations with increased confidence. Several benchmark ex- 
periments in differentially heated cavities have been conducted, such as Betts and Bokhari (2000) 
and Cooper et al. (2012). 


Advanced laser-based techniques for obtaining flow field measurements include Particle Image 
Velocimetry (PIV) and Laser Doppler Anemometry (LDA). Though such non-invasive techniques 
seem well suited to the kind of sensitive flows typical of passive cooling systems, they require 
optical access. This can be impossible for opaque fluids, but also difficult for water, since the pipes 
in real systems are likely to be made of opaque materials like metal. Transparent materials require 
optical properties matched to the fluid for best results in complex geometries (Hassan, 2019), and 
will typically have very different thermal properties to those likely to be used in any real system, 
which can interfere with the observed system performance, so need to be applied with care. 


In addition, the associated equipment can be cumbersome. For these reasons, it may be useful 
to use CFD modelling to identify areas where flow features of particular interest may occur and 
then focus the use of these techniques in these areas. CFD may also be used to indicate where 
additional measurements may be of value, such as using phased measurements in scenarios 
where the flow in the system is expected to be unstable. 


Velocity Measurement Techniques 


Pitot Tube: A pitot tube consists of two openings; one perpendicular to the flow, sampling local 
static pressure and one normal to the flow, sampling local total pressure (where the kinetic en- 
ergy of the flow is converted to a pressure increase). The difference in measured static and total 
pressure is used to calculate the local fluid flow velocity (Morris and Langari, 2020). 


Pitot tubes are a simple, cheap and widely used method for measuring velocity; however they have 
a low measurement accuracy. In addition, a pitot tube provides a single 1D velocity, although pitot 
cylinders (3-hole) and pitot spheres (5-hole) can measure 2D and 3D velocities respectively, but 
require calibration and may be bulky. 
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Constant Temperature Anemometry: Constant Temperature Anemometry (CTA), also known 
as Thermal Anemometry, is a technique for the measurement of 1D, 2D or 3D velocity and turbu- 
lence in gas and liquid flows, using hot-wire or hot-film probes inserted in the flow. CTA is partic- 
ularly suitable for the measurement of flows with very fast fluctuations at a point (high turbulence) 
and the study of flow micro-structures, where there is a need to resolve small flow eddies down to 
the order of tenths of a millimetre. 


Hot-wire anemometers consist of a thin electrically heated wire and the measurement principle is 
based on the cooling effect of a flow on a heated body, which reduces its resistance. Measurement 
of this resistance change is used to calculate the fluid velocity. Hot-wire anemometers have a fast 
response time (20 - 50 kHz); however they are not very robust due to the small diameter of the wire 
(Morris and Langari, 2020). Multiple wires are required to measure velocity components, which can 
become bulky and difficult to work with, so a single wire and traverse gear with phase-locking is 
often a more practical approach. 


Application areas include temperature, shear stress, velocity and turbulence measurements in e.g. 
jets, boundary layers and transitional flows. Wire sensors are used in gases and non-conducting 
liquids, while film sensors are primarily designed for use in water and other conducting liquids, as 
well as transition in gases. However, the temperature range for CTAs is limited, and they require 
careful, regular calibration. 


The Thermal Anemometry Grid Sensor is based on the same principle; however it contains mul- 
tiple temperature resistant resistors arranged in a grid. This enables spatially resolved velocity 
measurements to be obtained. 


Wire Mesh Sensor: Wire mesh sensors (WMSs) are an intrusive technique that allow the inves- 
tigation of multiphase flows with high spatial and temporal resolution. They can be used to obtain 
information about fluid velocity, void fraction, droplet/oubble size and distribution, interfacial area, 
film thickness, thermal distribution and flow regimes (Velasco Pefa and Rodriguez, 2015). 


Typically, a WMS consists of two parallel planes of wires; one plane containing transmitter wires, 
and one containing receiver wires. The planes are configured such that the wires on the top and 
bottom planes cross at an angle of 90°, forming a mesh grid of electrodes. This is placed into 
the cross-sectional area of the pipe or flow region of interest. The transmitter wires are activated 
sequentially, while the receiver wires are sampled in parallel. An electrical property (conductivity or 
permittivity) at each crossing point is evaluated to determine the fluid distribution across the cross- 
section. Velocity measurements are typically obtained using two WMS in different locations, and 
high accuracy 3D measurements can be made using two sensors separated by a few centimetres 
(Velasco Pefa and Rodriguez, 2015). 


WMS can be used at typical reactor temperatures and pressures, but is intrusive and so can influ- 
ence the flow field and size/shape of the bubbles. In addition, two phases are required in order to 
evaluate the conductivity/permittivity variation at each crossing e.g. air-water, steam-water or water 
with varying dissolved solute concentrations. This limits its applicability to velocity measurements, 
and therefore it is often used for measurements of mixing and void fraction. 
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Particle Image Velocimetry: PIV is a non-intrusive optical measurement technique which is 
used to obtain 2D or 3D fluid velocity measurements. PIV involves introducing small light scattering 
particles into the fluid flow. These particles are then illuminated by two consecutive pulses of a 
high intensity laser, and the scattered light across a plane of flow is recorded by a high resolution 
camera. The fluid velocity is then determined by analysis of the displacement of the particles 
between the successive frames. 


Stereo PIV measures three velocity components in a plane using two cameras, with the second 
camera at a different orientation such that it can record the third dimension component. Tomo- 
graphic PIV uses multiple cameras and can measure three velocity components in a volume. Time 
resolved PIV benefits from the advances in camera technology to acquire high resolution PIV im- 
ages at tens of thousands of frames per second with full camera resolution. 


In air flows, the seeding particles are typically oil drops in the range 1m to 5um. For water 
applications, the seeding is typically polystyrene, polyamide, hollow glass spheres or pepper in 
the range 5um to 100 um. It is worth noting that the seeding particles can affect the bulk fluid 
properties (Section 2.3), so it is important that the properties of the seeded fluid are understood 
and reported (and using in scaling assessments where appropriate, Section 3.3.1). 


Since PIV uses scattered light, the test fluid and test section walls must be optically transparent, 
which makes it difficult to use at reactor pressures and temperatures. PIV can be applied to flow 
regimes ranging from laminar, transitional to fully turbulent flows (Scarano, 2012). 


Laser Doppler Velocimetry/Anemometry: Laser Doppler Velocimetry (LDV) or LDA is a non- 
intrusive, optical technique for 1D, 2D and 3D point measurement of velocity and turbulence distri- 
bution in both free flows and internal flows. LDV requires that light scattering particles are present 
within the flowing fluid. The flow is then illuminated by a known frequency of laser, and the light 
scattered by the moving particles is detected by a photomultiplier. The difference in frequency be- 
tween the incident light and light reflected back is measured (i.e. the Doppler frequency), and used 
to calculate the velocity. 


LDV can be used in hostile environments and confined areas and can measure a wide range 
of velocities (Morris and Langari, 2020) with high spatial and temporal resolution. However, the 
method tends to be expensive and it can be difficult to collect data in close proximity to walls. 
Doppler Global Velocimetry is based on the same principles as LDV, but allows measurement of 
the velocity distribution within a plane. Since LDV uses lasers, the test fluid and test section walls 
must be optically transparent, which makes it difficult to use at reactor pressures and temperatures. 


Ultrasound Doppler Velocimetry: Ultrasound Doppler Velocimetry (UDV) is based on the pulse- 
echo method with velocity derived from shifts in positions between pulses, rather than shifts in 
frequency due to the Doppler effect. This requires acoustic inhomogeneities in the fluid to measure, 
which may be of natural origin as for many metal melts or artificial scattering particles have to be 
added. This technique is especially useful for opaque fluids or systems without optical access, 
since established optical flow measurement techniques such as PIV and LDV cannot be used. 


2D or 3D UDV can be used to measure two or three velocity components simultaneously along a 
line. UDV-2D is based on a three transducer system with one transducer as an emitter and two oth- 
ers as receivers. For UDV-3D measurements, a similar arrangement is used with three receivers. In 


67 of 91 


3.3.4 


3.3.5 


Methodology 


addition, spatial flow measurements can be achieved using several transducers arranged laterally 
to each other. However, the spatial and temporal resolution is limited by the number of transducers 
and sequential excitation required. 


Geometry, Boundary and Operating Conditions 


Clear and complete geometric information should be recorded, ideally with an ‘as-built’ Computer 
Aided Design (CAD) geometry. Where work is reported in a paper, geometric drawings dimen- 
sioned with sufficient detail to enable an accurate CAD representation to be extracted are important 
to enable following researchers to recreate the experiment. 


Seemingly minor geometric details that are often overlooked include bend radii and wall thick- 
nesses, which can vary due to flanges and manufacturing methods (e.g. bending pipes to form 
bends/coils can lead to thinner walls on the outer edge). The former is important since bends con- 
tribute to pressure losses within the pipe and inaccurate capture of these can lead to differences in 
mass flow rate and system performance at the low driving forces produced by natural convection. 
Wall thicknesses are important because the ability of the walls to absorb or release heat can affect 
the stability of the flow. For passive cooling systems, small differences can cause systems to switch 
between stable and unstable flow. 


It is equally important to obtain complete information about the boundary and operating conditions. 
For a passive cooling system, this would include at a minimum, the nature of the heat exchanger 
(i.e. best represented as a uniform heat flux or uniform temperature), the heat flux (or power) and, 
if appropriate, flow performance information for the secondary side, such as the mass flow rate 
and temperature. Values of ambient temperature should also be recorded and the details of any 
insulation surrounding any apparatus should be reported in detail, since this is often a source of 
uncertainty. In cases where the transient behaviour is of interest, detailed knowledge of any heater 
control systems is also important to record and report for future reference. 


Post-Processing 


There is often more than one way to compute key non-dimensional parameters that are used to 
characterise the system and, as such, the procedure used to compute the parameter should be 
recorded clearly, including quantifying the measurement uncertainty. The same procedure should 
be used for both numerical and experimental studies to ensure comparisons are valid (as discussed 
in Volume 1, Section 4.6.1). In addition, the temperature at which material properties have been 
evaluated at, and where or how that temperature was evaluated should also be reported. 


Material properties (Section 2.3) are often raised to various powers in non-dimensional parame- 
ters, meaning small differences in properties associated with small changes in temperatures can 
be greatly amplified and have a significant impact on the results. As an example, in a passive cool- 
ing system under natural circulation, the steady-state Re is often computed using a steady-state 
mass flow rate. Many experimental studies use an energy balance to compute this mass flow rate, 
since direct measurement is difficult. This can create uncertainty when comparing the Re from 
experimental work to values from numerical analysis. 
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To conclude the present volume, this section considers how the aspects presented on the preced- 
ing pages may change in the future. 


Over the past two decades, general purpose CFD codes (and particularly the RANS and URANS 
approaches contained within them) have reached a degree of maturity for the prediction of single- 
phase buoyant flows, and their use is now as routine as the longer established and very mature 
system codes. This maturity in development has been facilitated by a combination of extensive 
industrial use and modelling developments based on simple test cases. It is also expected that the 
use of coupled CFD and system code analyses will increase to better understand and predict the 
characteristics of areas of complex flow in parts of systems (Volume 2, Section 4). 


Looking further forward, Machine Learning (ML) could be used to drive analysis software in support 
of the more routine aspects of analysis, for example making decisions about how best to run 
previously developed and understood system and CFD models in support of a design sequence, 
or using a digital twin to predict the future performance of operating plant. In the nuclear industry 
however, the judgement of skilled engineers in understanding and justifying the safety of NPPs, to 
their own organisations, regulators and the public, is likely to remain of critical importance. 


In the short term, research and development is required on natural circulation, hybrid methods 
and turbulent heat transfer models to support the design of passive cooling systems in advanced 
reactors. 


Natural Circulation 


New reactor designs are making more extensive use of passive cooling systems for both normal 
operation and fault scenarios. Therefore, it is essential to be able to confidently predict the flow rate 
within the system loop under natural circulation conditions, and ensure that there is no chance of 
flow reversal or instability and the time taken for the flow to develop is understood. This challenge is 
highlighted in Section 2.4 and the current limitations of system codes is discussed in Section 3.1.1. 


Although CFD codes are able to more accurately predict the flow rate and onset of instability in 
natural circulation loops than system codes, the timescales for the flow to develop and extent of 
passive cooling systems means that using CFD analysis to support the design of reactor scale 
passive cooling systems is likely to be unfeasible. Therefore, further development is required to 
understand and develop best practice in this area. This should include: 


* A review of existing experimental data for natural circulation loops (Wilson, 2021) has been 
conducted as part of the Department for Business, Energy and Industrial Strategy (BEIS) 
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thermal hydraulics Nuclear Innovation Programme (NIP) (Volume 1, Appendix A). This demon- 
strated that there are very few high-fidelity measurements of natural circulation flows in repre- 
sentative reactor system loops. Therefore, there is a requirement for improved reactor scale 
natural circulation test facilities with high-fidelity measurements, such as non-intrusive opti- 
cal techniques providing instantaneous and time-averaged flow and thermal fields for model 
validation (Section 3.3.2, Section 3.3.3 and Volume 2, Section 3.5.1). 


The capability, limitations and adequacy of system codes to predict natural circulation in 
system loops needs to be better quantified and understood. This will enable best practice 
guidance to be developed in order to determine when a system code is suitable, when it 
should be coupled with a CFD code and when it is necessary to use CFD analysis. 


Though considered relatively mature, RANS and URANS are not well tested in complex 
passive cooling cases, such as modelling full natural circulation loops. Good progress has 
been made in Wilson (2021) as part of the BEIS thermal hydraulics NIP (Volume 1, Appendix 
A) to assess the benefit of using CFD for these applications. However, further work is required 
to apply and validate it on representative reactor-scale system loops. 


Hybrid methods 


As the cost of HPC resource reduces, so judgments on what is practicable for industrial modelling 
and the design process will also change. However, as outlined in Volumes 1, 2 and 3, the complexity 
and cost of moving to more detailed modelling approaches in CFD are significant. Therefore, while 
the use of LES is likely to increase slowly, it is unlikely to replace RANS (Hanjalic, 2005). 


The use of hybrid methods, such as WMLES, zonal LES, DES and SAS is also likely to in- 
crease. While significantly more computationally expensive than RANS, hybrid methods may offer 
increased fidelity and improved ability to resolve the flows in complex passive cooling systems, 
without the likely much greater cost of LES. In particular, this provides an alternative way of pro- 
viding validation or benchmark data as part of a CFD design optimisation. A typical design/safety 
justification process might involve the following steps: 


An initial validation exercise involving hybrid methods and experimental data relevant to the 
application of interest. 


This would provide a benchmark for validating simpler RANS/URANS simulations. 


The design/application could then be optimised through multiple RANS simulations with for- 
mal uncertainty quantification. 


Finally, a number of hybrid models could be run of the final solution to confirm the re- 
sults/conclusions from the RANS modelling. 


However, there is a requirement for more high-quality experimental data for detailed validation and 
performance assessment of hybrid methods for passive cooling systems. Further, more industry 
application is required in order to better understand the benefits and limitations of the current hybrid 
methods available (i.e. which method is most suitable for a particular application), and develop clear 
guidelines for appropriate mesh generation, model set-up and solution strategy. 
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Turbulent Heat Transfer Models 


Turbulent heat transfer models (Sections 2.2.4.4 and 3.2.6.5) simulate the turbulent heat trans- 
port, while RANS turbulence models simulate the turbulent momentum transport. These models 
have the potential to significantly improve the prediction of temperature gradients and turbulence 
generation for buoyancy affected flows (i.e. natural and mixed convection) using RANS turbulence 
models. 


This is an active area of research, particularly for low-Prandtl number flows of liquid metal (Volume 
5), although the options available in most CFD tools are limited. Examples include the implemen- 
tation of the SGDH and GGDH models in Fluent using User Defined Functions (UDFs) (Kumar 
and Dewan, 2013), and development of an explicit AHFM model in OpenFOAM (Manservisi and 
Menghini, 2014). Much of the recent development has been undertaken on implicit AHFMs within 
the Simulations and Experiments for the Safety Assessment of MEtal cooled reactors (SESAME) 
project. This includes the AHFM-2005 version that has been implemented in STAR-CCM+, and 
Nuclear Research and Consultancy Group (NRG)’s new variant of the model called AHFM-NRG+ 
(Shams et al., 2019). 


It is expected that these models will become more widely available within standard CFD analysis 
tools over the next few years (Volume 1, Section 4.5.4). Further validation will then be required to 
understand and assess their potential benefits for buoyancy driven, mixed and natural convection 
flows. 
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6 Nomenclature 


Latin Symbols 
A Area, m? 
At Atwood number (At = (p1 — p2)/(e1 + p2)) 
Bi Biot number (Bi = hL/kg) 
Cp, Cy Specific heat at constant pressure or volume, J kg~! K~* 
dorD_ Diameter (D, = 4A;;/Pcs for hydraulic diameter), m 
f Darcy-Weisbach friction factor 
Fo Fourier number (Fo = a,t/L?) 
Gr Grashof number (Gr = gL3Ap/v?p = gL?BAT/v?, using the Boussinesq approxi- 
mation Ap/p ~ —BAT, where AT is often taken as Ty, — Ts,00) 
g Acceleration due to gravity, ms~? 
h_ Specific enthalpy, J kg~+, Heat Transfer Coefficient (HTC), Wm~? K—! or height, m 
I Radiative intensity, W m~? sr~! or Wm~? sr~+ wm? for a spectral density, where sr 
(steradian) is solid angle 
J Radiosity, W m~? 
k Thermal conductivity, Wm-! K-? 
L Length or wall thickness, m 
M_ Molar mass of a species, kg kmol~? 
Ma _ Mach number (Ma = U/a, where a is the speed of sound) 
n_ Refractive index 
Nu Nusselt Number (Nu = AL/k¢) 
p Perimeter, m 
P Pressure (P; = static pressure, Py = total pressure), Nm~? or Pa 
Pe  Péclet number (Pe = RePr) 
Pr Prandtl number (Pr = cpu/ ke) 
q_ Heat flux (rate of heat transfer per unit area, g = Q/A), Wm~? 
Q_ Rate of heat transfer, W 
r Radius, m 
R_ Gas constant (for a particular gas, R = R/M), Jkg~!K~! 
R Universal gas constant, 8314.5 J kmol~+ K~+ 
Rip Thermal resistance, K W~! 
Ra Rayleigh number (Ra = GrPr) 
Re Reynolds number (Re = pUL/,, or for an internal flow Re = WD,,/A.<h) 
Ri Richardson number (Ri = Gr/Re’) 
Sr Strouhal number (Sr = fL/U, where f is frequency) 
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Nomenclature 'S, 


«x j<< 


Volume 3 


Stanton number (St = Nu/RePr) 

Time, s 

Temperature (7, = static temperature, 77 = total temperature), K 
Wall friction velocity (u, = \/Tw/p), ms~+ 

Velocity, ms~? or thermal transmittance, Wm? K~! 
Specific volume, m* kg~? 

Volume, m? 

Mass flow rate, kg s~+ 
Wall distance, m 


Non-dimensional wall distance (y* = yu, /v) 


Greek Symbols 


Qn SS FE YM KF Hm R WR 


4 


oa 


Thermal diffusivity (a = k/pcp), m?s~? 


Volumetric thermal expansion coefficient (8 = —(1/)(0p/0T)), K~* 
Ratio of specific heats (y = cp/c,) 

Emissivity or surface roughness height, m 

Absorption coefficient, m~? 

Wavelength, m 


Viscosity, kg m~1s~? 


Kinematic viscosity and momentum diffusivity (v = w/p), m?s~? 
Density, kg m~? 

Stefan Boltzmann constant, 5.67 x 10-8 W m-? K~* 

Shear stress, Nm~? 


Porosity or void fraction 


Subscripts and Modifications 


b 
cs 
f 


i 


wn 


8 


~ 


Bulk (mass-averaged) quantity 
Cross-sectional quantity 

Quantity relating to a fluid 

Quantity relating to a particular species 
Total (stagnation) quantity 

Turbulent quantity 

Static quantity or quantity relating to a solid 
Quantity relating to a wall or surface 
Quantity far from a wall or in free-stream 
Average quantity 

Molar quantity 

Varying quantity 
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7 Abbreviations 


AHFM 
ANT 
BEIS 
BEPU 
CAD 
CFD 
CFL 
CHF 
CHT 
CRP 
CSNI 
CTA 
DDES 
DES 
DFM 
DNS 
EBRSM 
ELES 
EMDAP 
EVM 
GEKO 
GGDH 
H2TS 
HPC 
HTC 
HTGR 
HVAG 
IAEA 
IAPWS 
IET 
ISP 
LCO 
LDA 
LDV 
LES 
LRR 
LWR 
ML 
NEA 
NIP 


Algebraic Heat Flux Model 

Advanced Nuclear Technology 

Department for Business, Energy and Industrial Strategy 
Best Estimate Plus Uncertainty 

Computer Aided Design 

Computational Fluid Dynamics 
Courant-Friedrichs Lewy 

Critical Heat Flux 

Conjugate Heat Transfer 

Coordinated Research Project 

Committee on the Safety of Nuclear Installations 
Constant Temperature Anemometry 

Delayed Detached Eddy Simulation 

Detached Eddy Simulation 

Differential Flux Model 

Direct Numerical Simulation 

Elliptic-Blending RSM 

Embedded LES 

Evaluation Methodology Development and Application Process 
Eddy Viscosity Model 

Generalized k-w 

Generalized Gradient Diffusion Hypothesis 
Hierarchical 2-Tiered Scaling 

High Performance Computing 

Heat Transfer Coefficient 

High Temperature Gas-cooled Reactor 

Heating, Ventilation, and Air Conditioning 
International Atomic Energy Agency 
International Association for the Properties of Water and Steam 
Integral Effect Test 

International Standard Problem 

Limits and Conditions for Operation 

Laser Doppler Anemometry 

Laser Doppler Velocimetry 

Large Eddy Simulation 

Launder-Reece-Rodi 

Light Water Reactor 

Machine Learning 

Nuclear Energy Agency 

Nuclear Innovation Programme 
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NPP Nuclear Power Plant 

NRG Nuclear Research and Consultancy Group 

NTH Nuclear Thermal Hydraulics 

O&M Operation and Maintenance 

OECD Organisation for Economic Co-operation and Development 

PIRT Phenomena Identification and Ranking Table 

PIV Particle Image Velocimetry 

PTS Pressurised Thermal Shock 

PWR Pressurised Water Reactor 

RANS Reynolds-Averaged Navier-Stokes 

RNG Re-Normalisation Group 

RPV Reactor Pressure Vessel 

RSM Reynolds Stress Model 

RVACS Reactor Vessel Auxiliary Cooling System 

SA Sensitivity Analysis 

SAS Scale-Adaptive Simulation 

SESAME __ Simulations and Experiments for the Safety Assessment of MEtal cooled reac- 

tors 

SGDH Simple Gradient Diffusion Hypothesis 

SGS Sub-Grid-Scale 

SSG Speziale-Sarkar-Gatski 

SST Shear Stress Transport 

TCL Two-Component Limit 

UDF User Defined Function 

UDV Ultrasound Doppler Velocimetry 

UQ Uncertainty Quantification 

URANS Unsteady Reynolds-Averaged Navier-Stokes 

US NRC United States Nuclear Regulatory Commission 

V&V Verification and Validation 

VVUQ Verification, Validation and Uncertainty Quantification 

WALE Wall-Adapting Local Eddy-Viscosity 

WMLES _ Wall Modeled Large Eddy Simulation 

WMS Wire mesh sensor 
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